source: trunk/Cbc/src/CbcBranchActual.cpp @ 912

Last change on this file since 912 was 912, checked in by ladanyi, 13 years ago

Incorporated changes from branches/heur

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 108.7 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 <cstdlib>
9#include <cmath>
10#include <cfloat>
11//#define CBC_DEBUG
12
13#include "CoinTypes.hpp"
14#include "OsiSolverInterface.hpp"
15#include "OsiSolverBranch.hpp"
16#include "CbcModel.hpp"
17#include "CbcMessage.hpp"
18#include "CbcBranchActual.hpp"
19#include "CoinSort.hpp"
20#include "CoinError.hpp"
21
22//##############################################################################
23
24// Default Constructor
25CbcClique::CbcClique ()
26  : CbcObject(),
27    numberMembers_(0),
28    numberNonSOSMembers_(0),
29    members_(NULL),
30    type_(NULL),
31    cliqueType_(-1),
32    slack_(-1)
33{
34}
35
36// Useful constructor (which are integer indices)
37CbcClique::CbcClique (CbcModel * model, int cliqueType, int numberMembers,
38           const int * which, const char * type, int identifier,int slack)
39  : CbcObject(model)
40{
41  id_=identifier;
42  numberMembers_=numberMembers;
43  if (numberMembers_) {
44    members_ = new int[numberMembers_];
45    memcpy(members_,which,numberMembers_*sizeof(int));
46    type_ = new char[numberMembers_];
47    if (type) {
48      memcpy(type_,type,numberMembers_*sizeof(char));
49    } else {
50      for (int i=0;i<numberMembers_;i++)
51        type_[i]=1;
52    }
53  } else {
54    members_ = NULL;
55    type_ = NULL;
56  }
57  // Find out how many non sos
58  int i;
59  numberNonSOSMembers_=0;
60  for (i=0;i<numberMembers_;i++)
61    if (!type_[i])
62      numberNonSOSMembers_++;
63  cliqueType_ = cliqueType;
64  slack_ = slack;
65}
66
67// Copy constructor
68CbcClique::CbcClique ( const CbcClique & rhs)
69  :CbcObject(rhs)
70{
71  numberMembers_ = rhs.numberMembers_;
72  numberNonSOSMembers_ = rhs.numberNonSOSMembers_;
73  if (numberMembers_) {
74    members_ = new int[numberMembers_];
75    memcpy(members_,rhs.members_,numberMembers_*sizeof(int));
76    type_ = new char[numberMembers_];
77    memcpy(type_,rhs.type_,numberMembers_*sizeof(char));
78  } else {
79    members_ = NULL;
80    type_ = NULL;
81  }
82  cliqueType_ = rhs.cliqueType_;
83  slack_ = rhs.slack_;
84}
85
86// Clone
87CbcObject *
88CbcClique::clone() const
89{
90  return new CbcClique(*this);
91}
92
93// Assignment operator
94CbcClique & 
95CbcClique::operator=( const CbcClique& rhs)
96{
97  if (this!=&rhs) {
98    CbcObject::operator=(rhs);
99    delete [] members_;
100    delete [] type_;
101    numberMembers_ = rhs.numberMembers_;
102    numberNonSOSMembers_ = rhs.numberNonSOSMembers_;
103    if (numberMembers_) {
104      members_ = new int[numberMembers_];
105      memcpy(members_,rhs.members_,numberMembers_*sizeof(int));
106      type_ = new char[numberMembers_];
107      memcpy(type_,rhs.type_,numberMembers_*sizeof(char));
108    } else {
109      members_ = NULL;
110      type_ = NULL;
111    }
112    cliqueType_ = rhs.cliqueType_;
113    slack_ = rhs.slack_;
114  }
115  return *this;
116}
117
118// Destructor
119CbcClique::~CbcClique ()
120{
121  delete [] members_;
122  delete [] type_;
123}
124
125// Infeasibility - large is 0.5
126double 
127CbcClique::infeasibility(int & preferredWay) const
128{
129  int numberUnsatis=0, numberFree=0;
130  int j;
131  const int * integer = model_->integerVariable();
132  OsiSolverInterface * solver = model_->solver();
133  const double * solution = model_->testSolution();
134  const double * lower = solver->getColLower();
135  const double * upper = solver->getColUpper();
136  double largestValue=0.0;
137  double integerTolerance = 
138    model_->getDblParam(CbcModel::CbcIntegerTolerance);
139  double * sort = new double[numberMembers_];
140
141  double slackValue=0.0;
142  for (j=0;j<numberMembers_;j++) {
143    int sequence = members_[j];
144    int iColumn = integer[sequence];
145    double value = solution[iColumn];
146    value = CoinMax(value, lower[iColumn]);
147    value = CoinMin(value, upper[iColumn]);
148    double nearest = floor(value+0.5);
149    double distance = fabs(value-nearest);
150    if (distance>integerTolerance) {
151      if (!type_[j])
152        value = 1.0-value; // non SOS
153      // if slack then choose that
154      if (j==slack_&&value>0.05)
155        slackValue = value;
156      largestValue = CoinMax(value,largestValue);
157      sort[numberUnsatis++]=-value;
158    } else if (upper[iColumn]>lower[iColumn]) {
159      numberFree++;
160    }
161  }
162  preferredWay=1;
163  double otherWay = 0.0;
164  if (numberUnsatis) {
165    // sort
166    std::sort(sort,sort+numberUnsatis);
167    for (j=0;j<numberUnsatis;j++) {
168      if ((j&1)!=0)
169        otherWay += -sort[j];
170    }
171    // Need to think more
172    double value = 0.2*numberUnsatis+0.01*(numberMembers_-numberFree);
173    if (fabs(largestValue-0.5)<0.1) {
174      // close to half
175      value +=0.1;
176    }
177    if (slackValue) {
178      // branching on slack
179      value += slackValue;
180    }
181    // scale other way
182    otherWay *= value/(1.0-otherWay);
183    delete [] sort;
184    return value;
185  } else {
186    delete [] sort;
187    return 0.0; // satisfied
188  }
189}
190
191// This looks at solution and sets bounds to contain solution
192void 
193CbcClique::feasibleRegion()
194{
195  int j;
196  const int * integer = model_->integerVariable();
197  OsiSolverInterface * solver = model_->solver();
198  const double * solution = model_->testSolution();
199  const double * lower = solver->getColLower();
200  const double * upper = solver->getColUpper();
201#ifndef NDEBUG
202  double integerTolerance = 
203    model_->getDblParam(CbcModel::CbcIntegerTolerance);
204#endif 
205  for (j=0;j<numberMembers_;j++) {
206    int sequence = members_[j];
207    int iColumn = integer[sequence];
208    double value = solution[iColumn];
209    value = CoinMax(value, lower[iColumn]);
210    value = CoinMin(value, upper[iColumn]);
211    double nearest = floor(value+0.5);
212#ifndef NDEBUG
213    double distance = fabs(value-nearest);
214    assert(distance<=integerTolerance);
215#endif
216    solver->setColLower(iColumn,nearest);
217    solver->setColUpper(iColumn,nearest);
218  }
219}
220// Redoes data when sequence numbers change
221void 
222CbcClique::redoSequenceEtc(CbcModel * model, int numberColumns, const int * originalColumns)
223{
224  model_=model;
225  int n2=0;
226  for (int j=0;j<numberMembers_;j++) {
227    int iColumn = members_[j];
228    int i;
229    for (i=0;i<numberColumns;i++) {
230      if (originalColumns[i]==iColumn)
231        break;
232    }
233    if (i<numberColumns) {
234      members_[n2]=i;
235      type_[n2++]=type_[j];
236    }
237  }
238  if (n2<numberMembers_) {
239    //printf("** SOS number of members reduced from %d to %d!\n",numberMembers_,n2);
240    numberMembers_=n2;
241  }
242  // Find out how many non sos
243  int i;
244  numberNonSOSMembers_=0;
245  for (i=0;i<numberMembers_;i++)
246    if (!type_[i])
247      numberNonSOSMembers_++;
248}
249
250
251// Creates a branching object
252CbcBranchingObject * 
253CbcClique::createBranch(int way) 
254{
255  int numberUnsatis=0;
256  int j;
257  int nUp=0;
258  int nDown=0;
259  int numberFree=numberMembers_;
260  const int * integer = model_->integerVariable();
261  OsiSolverInterface * solver = model_->solver();
262  const double * solution = model_->testSolution();
263  const double * lower = solver->getColLower();
264  const double * upper = solver->getColUpper();
265  int * upList = new int[numberMembers_];
266  int * downList = new int[numberMembers_];
267  double * sort = new double[numberMembers_];
268  double integerTolerance = 
269      model_->getDblParam(CbcModel::CbcIntegerTolerance);
270
271  double slackValue=0.0;
272  for (j=0;j<numberMembers_;j++) {
273    int sequence = members_[j];
274    int iColumn = integer[sequence];
275    double value = solution[iColumn];
276    value = CoinMax(value, lower[iColumn]);
277    value = CoinMin(value, upper[iColumn]);
278    double nearest = floor(value+0.5);
279    double distance = fabs(value-nearest);
280    if (distance>integerTolerance) {
281      if (!type_[j])
282        value = 1.0-value; // non SOS
283      // if slack then choose that
284      if (j==slack_&&value>0.05)
285        slackValue = value;
286      value = -value; // for sort
287      upList[numberUnsatis]=j;
288      sort[numberUnsatis++]=value;
289    } else if (upper[iColumn]>lower[iColumn]) {
290      upList[--numberFree]=j;
291    }
292  }
293  assert (numberUnsatis);
294  if (!slackValue) {
295    // sort
296    CoinSort_2(sort,sort+numberUnsatis,upList);
297    // put first in up etc
298    int kWay=1;
299    for (j=0;j<numberUnsatis;j++) {
300      if (kWay>0) 
301        upList[nUp++]=upList[j];
302      else
303        downList[nDown++]=upList[j];
304      kWay = -kWay;
305    }
306    for (j=numberFree;j<numberMembers_;j++) {
307      if (kWay>0)
308        upList[nUp++]=upList[j];
309      else
310        downList[nDown++]=upList[j];
311      kWay = -kWay;
312    }
313  } else {
314    // put slack to 0 in first way
315    nUp = 1;
316    upList[0]=slack_;
317    for (j=0;j<numberUnsatis;j++) {
318      downList[nDown++]=upList[j];
319    }
320    for (j=numberFree;j<numberMembers_;j++) {
321      downList[nDown++]=upList[j];
322    }
323  }
324  // create object
325  CbcBranchingObject * branch;
326  if (numberMembers_ <=64)
327     branch = new CbcCliqueBranchingObject(model_,this,way,
328                                         nDown,downList,nUp,upList);
329  else
330    branch = new CbcLongCliqueBranchingObject(model_,this,way,
331                                            nDown,downList,nUp,upList);
332  delete [] upList;
333  delete [] downList;
334  delete [] sort;
335  return branch;
336}
337
338//##############################################################################
339
340// Default Constructor
341CbcSOS::CbcSOS ()
342  : CbcObject(),
343    members_(NULL),
344    weights_(NULL),
345    numberMembers_(0),
346    sosType_(-1),
347    integerValued_(false)
348{
349}
350
351// Useful constructor (which are indices)
352CbcSOS::CbcSOS (CbcModel * model,  int numberMembers,
353           const int * which, const double * weights, int identifier,int type)
354  : CbcObject(model),
355    numberMembers_(numberMembers),
356    sosType_(type)
357{
358  id_=identifier;
359  integerValued_ = type==1;
360  if (integerValued_) {
361    // check all members integer
362    OsiSolverInterface * solver = model->solver();
363    if (solver) {
364      for (int i=0;i<numberMembers_;i++) {
365        if (!solver->isInteger(which[i]))
366          integerValued_=false;
367      }
368    } else {
369      // can't tell
370      integerValued_=false;
371    }
372  }
373  if (numberMembers_) {
374    members_ = new int[numberMembers_];
375    weights_ = new double[numberMembers_];
376    memcpy(members_,which,numberMembers_*sizeof(int));
377    if (weights) {
378      memcpy(weights_,weights,numberMembers_*sizeof(double));
379    } else {
380      for (int i=0;i<numberMembers_;i++)
381        weights_[i]=i;
382    }
383    // sort so weights increasing
384    CoinSort_2(weights_,weights_+numberMembers_,members_);
385    double last = -COIN_DBL_MAX;
386    int i;
387    for (i=0;i<numberMembers_;i++) {
388      double possible = CoinMax(last+1.0e-10,weights_[i]);
389      weights_[i] = possible;
390      last=possible;
391    }
392  } else {
393    members_ = NULL;
394    weights_ = NULL;
395  }
396  assert (sosType_>0&&sosType_<3);
397}
398
399// Copy constructor
400CbcSOS::CbcSOS ( const CbcSOS & rhs)
401  :CbcObject(rhs)
402{
403  numberMembers_ = rhs.numberMembers_;
404  sosType_ = rhs.sosType_;
405  integerValued_ = rhs.integerValued_;
406  if (numberMembers_) {
407    members_ = new int[numberMembers_];
408    weights_ = new double[numberMembers_];
409    memcpy(members_,rhs.members_,numberMembers_*sizeof(int));
410    memcpy(weights_,rhs.weights_,numberMembers_*sizeof(double));
411  } else {
412    members_ = NULL;
413    weights_ = NULL;
414  }
415}
416
417// Clone
418CbcObject *
419CbcSOS::clone() const
420{
421  return new CbcSOS(*this);
422}
423
424// Assignment operator
425CbcSOS & 
426CbcSOS::operator=( const CbcSOS& rhs)
427{
428  if (this!=&rhs) {
429    CbcObject::operator=(rhs);
430    delete [] members_;
431    delete [] weights_;
432    numberMembers_ = rhs.numberMembers_;
433    sosType_ = rhs.sosType_;
434    integerValued_ = rhs.integerValued_;
435    if (numberMembers_) {
436      members_ = new int[numberMembers_];
437      weights_ = new double[numberMembers_];
438      memcpy(members_,rhs.members_,numberMembers_*sizeof(int));
439      memcpy(weights_,rhs.weights_,numberMembers_*sizeof(double));
440    } else {
441      members_ = NULL;
442      weights_ = NULL;
443    }
444  }
445  return *this;
446}
447
448// Destructor
449CbcSOS::~CbcSOS ()
450{
451  delete [] members_;
452  delete [] weights_;
453}
454
455// Infeasibility - large is 0.5
456double 
457CbcSOS::infeasibility(int & preferredWay) const
458{
459  int j;
460  int firstNonZero=-1;
461  int lastNonZero = -1;
462  OsiSolverInterface * solver = model_->solver();
463  const double * solution = model_->testSolution();
464  //const double * lower = solver->getColLower();
465  const double * upper = solver->getColUpper();
466  //double largestValue=0.0;
467  double integerTolerance = 
468    model_->getDblParam(CbcModel::CbcIntegerTolerance);
469  double weight = 0.0;
470  double sum =0.0;
471
472  // check bounds etc
473  double lastWeight=-1.0e100;
474  for (j=0;j<numberMembers_;j++) {
475    int iColumn = members_[j];
476    if (lastWeight>=weights_[j]-1.0e-7)
477      throw CoinError("Weights too close together in SOS","infeasibility","CbcSOS");
478    double value = CoinMax(0.0,solution[iColumn]);
479    sum += value;
480    if (value>integerTolerance&&upper[iColumn]) {
481      // Possibly due to scaling a fixed variable might slip through
482      if (value>upper[iColumn]) {
483        value=upper[iColumn];
484        // Could change to #ifdef CBC_DEBUG
485#ifndef NDEBUG
486        if (model_->messageHandler()->logLevel()>2)
487          printf("** Variable %d (%d) has value %g and upper bound of %g\n",
488                 iColumn,j,value,upper[iColumn]);
489#endif
490      } 
491      weight += weights_[j]*value;
492      if (firstNonZero<0)
493        firstNonZero=j;
494      lastNonZero=j;
495    }
496  }
497  preferredWay=1;
498  if (lastNonZero-firstNonZero>=sosType_) {
499    // find where to branch
500    assert (sum>0.0);
501    weight /= sum;
502    //int iWhere;
503    //for (iWhere=firstNonZero;iWhere<lastNonZero;iWhere++)
504    //if (weight<weights_[iWhere+1])
505    //break;
506    // probably best to use pseudo duals
507    double value = lastNonZero-firstNonZero+1;
508    value *= 0.5/((double) numberMembers_);
509    return value;
510  } else {
511    return 0.0; // satisfied
512  }
513}
514
515// This looks at solution and sets bounds to contain solution
516void 
517CbcSOS::feasibleRegion()
518{
519  int j;
520  int firstNonZero=-1;
521  int lastNonZero = -1;
522  OsiSolverInterface * solver = model_->solver();
523  const double * solution = model_->testSolution();
524  //const double * lower = solver->getColLower();
525  const double * upper = solver->getColUpper();
526  double integerTolerance = 
527    model_->getDblParam(CbcModel::CbcIntegerTolerance);
528  double weight = 0.0;
529  double sum =0.0;
530
531  for (j=0;j<numberMembers_;j++) {
532    int iColumn = members_[j];
533    double value = CoinMax(0.0,solution[iColumn]);
534    sum += value;
535    if (value>integerTolerance&&upper[iColumn]) {
536      weight += weights_[j]*value;
537      if (firstNonZero<0)
538        firstNonZero=j;
539      lastNonZero=j;
540    }
541  }
542  assert (lastNonZero-firstNonZero<sosType_) ;
543  for (j=0;j<firstNonZero;j++) {
544    int iColumn = members_[j];
545    solver->setColUpper(iColumn,0.0);
546  }
547  for (j=lastNonZero+1;j<numberMembers_;j++) {
548    int iColumn = members_[j];
549    solver->setColUpper(iColumn,0.0);
550  }
551}
552// Redoes data when sequence numbers change
553void 
554CbcSOS::redoSequenceEtc(CbcModel * model, int numberColumns, const int * originalColumns)
555{
556  model_=model;
557  int n2=0;
558  for (int j=0;j<numberMembers_;j++) {
559    int iColumn = members_[j];
560    int i;
561    for (i=0;i<numberColumns;i++) {
562      if (originalColumns[i]==iColumn)
563        break;
564    }
565    if (i<numberColumns) {
566      members_[n2]=i;
567      weights_[n2++]=weights_[j];
568    }
569  }
570  if (n2<numberMembers_) {
571    //printf("** SOS number of members reduced from %d to %d!\n",numberMembers_,n2);
572    numberMembers_=n2;
573  }
574}
575
576
577// Creates a branching object
578CbcBranchingObject * 
579CbcSOS::createBranch(int way) 
580{
581  int j;
582  const double * solution = model_->testSolution();
583  double integerTolerance = 
584      model_->getDblParam(CbcModel::CbcIntegerTolerance);
585  OsiSolverInterface * solver = model_->solver();
586  const double * upper = solver->getColUpper();
587  int firstNonFixed=-1;
588  int lastNonFixed=-1;
589  int firstNonZero=-1;
590  int lastNonZero = -1;
591  double weight = 0.0;
592  double sum =0.0;
593  for (j=0;j<numberMembers_;j++) {
594    int iColumn = members_[j];
595    if (upper[iColumn]) {
596      double value = CoinMax(0.0,solution[iColumn]);
597      sum += value;
598      if (firstNonFixed<0)
599        firstNonFixed=j;
600      lastNonFixed=j;
601      if (value>integerTolerance) {
602        weight += weights_[j]*value;
603        if (firstNonZero<0)
604          firstNonZero=j;
605        lastNonZero=j;
606      }
607    }
608  }
609  assert (lastNonZero-firstNonZero>=sosType_) ;
610  // find where to branch
611  assert (sum>0.0);
612  weight /= sum;
613  int iWhere;
614  double separator=0.0;
615  for (iWhere=firstNonZero;iWhere<lastNonZero;iWhere++) 
616    if (weight<weights_[iWhere+1])
617      break;
618  if (sosType_==1) {
619    // SOS 1
620    separator = 0.5 *(weights_[iWhere]+weights_[iWhere+1]);
621  } else {
622    // SOS 2
623    if (iWhere==firstNonFixed)
624      iWhere++;;
625    if (iWhere==lastNonFixed-1)
626      iWhere = lastNonFixed-2;
627    separator = weights_[iWhere+1];
628  }
629  // create object
630  CbcBranchingObject * branch;
631  branch = new CbcSOSBranchingObject(model_,this,way,separator);
632  return branch;
633}
634
635/* Create an OsiSolverBranch object
636   
637This returns NULL if branch not represented by bound changes
638*/
639OsiSolverBranch * 
640CbcSOS::solverBranch() const
641{
642  int j;
643  const double * solution = model_->testSolution();
644  double integerTolerance = 
645      model_->getDblParam(CbcModel::CbcIntegerTolerance);
646  OsiSolverInterface * solver = model_->solver();
647  const double * upper = solver->getColUpper();
648  int firstNonFixed=-1;
649  int lastNonFixed=-1;
650  int firstNonZero=-1;
651  int lastNonZero = -1;
652  double weight = 0.0;
653  double sum =0.0;
654  double * fix = new double[numberMembers_];
655  int * which = new int[numberMembers_];
656  for (j=0;j<numberMembers_;j++) {
657    int iColumn = members_[j];
658    // fix all on one side or other (even if fixed)
659    fix[j]=0.0;
660    which[j]=iColumn;
661    if (upper[iColumn]) {
662      double value = CoinMax(0.0,solution[iColumn]);
663      sum += value;
664      if (firstNonFixed<0)
665        firstNonFixed=j;
666      lastNonFixed=j;
667      if (value>integerTolerance) {
668        weight += weights_[j]*value;
669        if (firstNonZero<0)
670          firstNonZero=j;
671        lastNonZero=j;
672      }
673    }
674  }
675  assert (lastNonZero-firstNonZero>=sosType_) ;
676  // find where to branch
677  assert (sum>0.0);
678  weight /= sum;
679  // down branch fixes ones above weight to 0
680  int iWhere;
681  int iDownStart=0;
682  int iUpEnd=0;
683  for (iWhere=firstNonZero;iWhere<lastNonZero;iWhere++) 
684    if (weight<weights_[iWhere+1])
685      break;
686  if (sosType_==1) {
687    // SOS 1
688    iUpEnd=iWhere+1;
689    iDownStart=iUpEnd;
690  } else {
691    // SOS 2
692    if (iWhere==firstNonFixed)
693      iWhere++;;
694    if (iWhere==lastNonFixed-1)
695      iWhere = lastNonFixed-2;
696    iUpEnd=iWhere+1;
697    iDownStart=iUpEnd+1;
698  }
699  //
700  OsiSolverBranch * branch = new OsiSolverBranch();
701  branch->addBranch(-1,0,NULL,NULL,numberMembers_-iDownStart,which+iDownStart,fix);
702  branch->addBranch(1,0,NULL,NULL,iUpEnd,which,fix);
703  delete [] fix;
704  delete [] which;
705  return branch;
706}
707// Construct an OsiSOS object
708OsiSOS * 
709CbcSOS::osiObject(const OsiSolverInterface * solver) const
710{
711  OsiSOS * obj = new OsiSOS(solver,numberMembers_,members_,weights_,sosType_);
712  obj->setPriority(priority());
713  return obj;
714}
715
716//##############################################################################
717
718/** Default Constructor
719
720  Equivalent to an unspecified binary variable.
721*/
722CbcSimpleInteger::CbcSimpleInteger ()
723  : CbcObject(),
724    originalLower_(0.0),
725    originalUpper_(1.0),
726    breakEven_(0.5),
727    columnNumber_(-1),
728    preferredWay_(0)
729{
730}
731
732/** Useful constructor
733
734  Loads actual upper & lower bounds for the specified variable.
735*/
736CbcSimpleInteger::CbcSimpleInteger ( CbcModel * model, int iColumn, double breakEven)
737  : CbcObject(model)
738{
739  columnNumber_ = iColumn ;
740  originalLower_ = model->solver()->getColLower()[columnNumber_] ;
741  originalUpper_ = model->solver()->getColUpper()[columnNumber_] ;
742  breakEven_ = breakEven;
743  assert (breakEven_>0.0&&breakEven_<1.0);
744  preferredWay_ = 0;
745}
746
747
748// Copy constructor
749CbcSimpleInteger::CbcSimpleInteger ( const CbcSimpleInteger & rhs)
750  :CbcObject(rhs)
751
752{
753  columnNumber_ = rhs.columnNumber_;
754  originalLower_ = rhs.originalLower_;
755  originalUpper_ = rhs.originalUpper_;
756  breakEven_ = rhs.breakEven_;
757  preferredWay_ = rhs.preferredWay_;
758}
759
760// Clone
761CbcObject *
762CbcSimpleInteger::clone() const
763{
764  return new CbcSimpleInteger(*this);
765}
766
767// Assignment operator
768CbcSimpleInteger & 
769CbcSimpleInteger::operator=( const CbcSimpleInteger& rhs)
770{
771  if (this!=&rhs) {
772    CbcObject::operator=(rhs);
773    columnNumber_ = rhs.columnNumber_;
774    originalLower_ = rhs.originalLower_;
775    originalUpper_ = rhs.originalUpper_;
776    breakEven_ = rhs.breakEven_;
777    preferredWay_ = rhs.preferredWay_;
778  }
779  return *this;
780}
781
782// Destructor
783CbcSimpleInteger::~CbcSimpleInteger ()
784{
785}
786// Construct an OsiSimpleInteger object
787OsiSimpleInteger * 
788CbcSimpleInteger::osiObject() const
789{
790  OsiSimpleInteger * obj = new OsiSimpleInteger(columnNumber_,
791                                                originalLower_,originalUpper_);
792  obj->setPriority(priority());
793  return obj;
794}
795
796double
797CbcSimpleInteger::infeasibility(const OsiSolverInterface * solver, const OsiBranchingInformation * info,
798                         int & preferredWay) const
799{
800  double value = info->solution_[columnNumber_];
801  value = CoinMax(value, info->lower_[columnNumber_]);
802  value = CoinMin(value, info->upper_[columnNumber_]);
803  double nearest = floor(value+(1.0-breakEven_));
804  assert (breakEven_>0.0&&breakEven_<1.0);
805  if (nearest>value) 
806    preferredWay=1;
807  else
808    preferredWay=-1;
809  if (preferredWay_)
810    preferredWay=preferredWay_;
811  double weight = fabs(value-nearest);
812  // normalize so weight is 0.5 at break even
813  if (nearest<value)
814    weight = (0.5/breakEven_)*weight;
815  else
816    weight = (0.5/(1.0-breakEven_))*weight;
817  if (fabs(value-nearest)<=info->integerTolerance_) 
818    return 0.0;
819  else
820    return weight;
821}
822double 
823CbcSimpleInteger::feasibleRegion(OsiSolverInterface * solver, const OsiBranchingInformation * info) const 
824{
825  double value = info->solution_[columnNumber_];
826  double newValue = CoinMax(value, info->lower_[columnNumber_]);
827  newValue = CoinMin(newValue, info->upper_[columnNumber_]);
828  newValue = floor(newValue+0.5);
829  solver->setColLower(columnNumber_,newValue);
830  solver->setColUpper(columnNumber_,newValue);
831  return fabs(value-newValue);
832}
833
834/* Create an OsiSolverBranch object
835   
836This returns NULL if branch not represented by bound changes
837*/
838OsiSolverBranch * 
839CbcSimpleInteger::solverBranch(OsiSolverInterface * solver, const OsiBranchingInformation * info) const
840{
841  double value = info->solution_[columnNumber_];
842  value = CoinMax(value, info->lower_[columnNumber_]);
843  value = CoinMin(value, info->upper_[columnNumber_]);
844  assert (info->upper_[columnNumber_]>info->lower_[columnNumber_]);
845#ifndef NDEBUG
846  double nearest = floor(value+0.5);
847  assert (fabs(value-nearest)>info->integerTolerance_);
848#endif
849  OsiSolverBranch * branch = new OsiSolverBranch();
850  branch->addBranch(columnNumber_,value);
851  return branch;
852}
853// Creates a branching object
854CbcBranchingObject * 
855CbcSimpleInteger::createBranch(OsiSolverInterface * solver, const OsiBranchingInformation * info, int way) 
856{
857  CbcIntegerBranchingObject * branch = new CbcIntegerBranchingObject(model_,0,-1,0.5);
858  fillCreateBranch(branch,info,way);
859  return branch;
860}
861// Fills in a created branching object
862void 
863CbcSimpleInteger::fillCreateBranch(CbcIntegerBranchingObject * branch, const OsiBranchingInformation * info, int way) 
864{
865  branch->setOriginalObject(this);
866  double value = info->solution_[columnNumber_];
867  value = CoinMax(value, info->lower_[columnNumber_]);
868  value = CoinMin(value, info->upper_[columnNumber_]);
869  assert (info->upper_[columnNumber_]>info->lower_[columnNumber_]);
870  if (!info->hotstartSolution_) {
871#ifndef NDEBUG
872    double nearest = floor(value+0.5);
873    assert (fabs(value-nearest)>info->integerTolerance_);
874#endif
875  } else {
876    double targetValue = info->hotstartSolution_[columnNumber_];
877    if (way>0)
878      value = targetValue-0.1;
879    else
880      value = targetValue+0.1;
881  }
882  branch->fillPart(columnNumber_,way,value);
883}
884/* Column number if single column object -1 otherwise,
885   so returns >= 0
886   Used by heuristics
887*/
888int 
889CbcSimpleInteger::columnNumber() const
890{
891  return columnNumber_;
892}
893/* Reset variable bounds to their original values.
894 
895    Bounds may be tightened, so it may be good to be able to set this info in object.
896*/
897void 
898CbcSimpleInteger::resetBounds(const OsiSolverInterface * solver) 
899{
900  originalLower_ = solver->getColLower()[columnNumber_] ;
901  originalUpper_ = solver->getColUpper()[columnNumber_] ;
902}
903
904/*  Change column numbers after preprocessing
905 */
906void 
907CbcSimpleInteger::resetSequenceEtc(int numberColumns, const int * originalColumns) 
908{
909  assert (numberColumns>0);
910  int iColumn;
911#if 0
912  for (iColumn=0;iColumn<numberColumns;iColumn++) {
913    if (columnNumber_==originalColumns[iColumn])
914      break;
915  }
916  assert (iColumn<numberColumns);
917#else
918  iColumn=originalColumns[columnNumber_];
919  assert (iColumn>=0);
920#endif
921  columnNumber_ = iColumn;
922}
923
924// Infeasibility - large is 0.5
925double 
926CbcSimpleInteger::infeasibility(int & preferredWay) const
927{
928  OsiBranchingInformation info(model_->solver(),model_->normalSolver(),false);
929  return infeasibility(model_->solver(),&info,preferredWay);
930}
931
932// This looks at solution and sets bounds to contain solution
933/** More precisely: it first forces the variable within the existing
934    bounds, and then tightens the bounds to fix the variable at the
935    nearest integer value.
936*/
937void 
938CbcSimpleInteger::feasibleRegion()
939{
940  abort();
941}
942CbcBranchingObject * 
943CbcSimpleInteger::createBranch( int way) 
944{
945  abort();
946  return NULL;
947}
948
949//##############################################################################
950
951// Default Constructor
952CbcIntegerBranchingObject::CbcIntegerBranchingObject()
953  :CbcBranchingObject()
954{
955  down_[0] = 0.0;
956  down_[1] = 0.0;
957  up_[0] = 0.0;
958  up_[1] = 0.0;
959#ifdef FUNNY_BRANCHING
960  variables_ = NULL;
961  newBounds_ = NULL;
962  numberExtraChangedBounds_ = 0;
963#endif
964}
965// Useful constructor
966CbcIntegerBranchingObject::CbcIntegerBranchingObject (CbcModel * model, 
967                                                      int variable, int way , double value)
968  :CbcBranchingObject(model,variable,way,value)
969{
970  int iColumn = variable;
971  assert (model_->solver()->getNumCols()>0);
972  down_[0] = model_->solver()->getColLower()[iColumn];
973  down_[1] = floor(value_);
974  up_[0] = ceil(value_);
975  up_[1] = model->getColUpper()[iColumn];
976#ifdef FUNNY_BRANCHING
977  variables_ = NULL;
978  newBounds_ = NULL;
979  numberExtraChangedBounds_ = 0;
980#endif
981}
982// Does part of constructor
983void 
984CbcIntegerBranchingObject::fillPart (int variable,
985                                 int way , double value) 
986{
987  //originalObject_=NULL;
988  branchIndex_=0;
989  value_=value;
990  numberBranches_=2;
991  //model_= model;
992  //originalCbcObject_=NULL;
993  variable_=variable;
994  way_=way;
995  int iColumn = variable;
996  down_[0] = model_->solver()->getColLower()[iColumn];
997  down_[1] = floor(value_);
998  up_[0] = ceil(value_);
999  up_[1] = model_->getColUpper()[iColumn];
1000}
1001// Useful constructor for fixing
1002CbcIntegerBranchingObject::CbcIntegerBranchingObject (CbcModel * model, 
1003                                                      int variable, int way,
1004                                                      double lowerValue, 
1005                                                      double upperValue)
1006  :CbcBranchingObject(model,variable,way,lowerValue)
1007{
1008  setNumberBranchesLeft(1);
1009  down_[0] = lowerValue;
1010  down_[1] = upperValue;
1011  up_[0] = lowerValue;
1012  up_[1] = upperValue;
1013#ifdef FUNNY_BRANCHING
1014  variables_ = NULL;
1015  newBounds_ = NULL;
1016  numberExtraChangedBounds_ = 0;
1017#endif
1018}
1019 
1020
1021// Copy constructor
1022CbcIntegerBranchingObject::CbcIntegerBranchingObject ( const CbcIntegerBranchingObject & rhs) :CbcBranchingObject(rhs)
1023{
1024  down_[0] = rhs.down_[0];
1025  down_[1] = rhs.down_[1];
1026  up_[0] = rhs.up_[0];
1027  up_[1] = rhs.up_[1];
1028#ifdef FUNNY_BRANCHING
1029  numberExtraChangedBounds_ = rhs.numberExtraChangedBounds_;
1030  int size = numberExtraChangedBounds_*(sizeof(double)+sizeof(int));
1031  char * temp = new char [size];
1032  newBounds_ = (double *) temp;
1033  variables_ = (int *) (newBounds_+numberExtraChangedBounds_);
1034
1035  int i ;
1036  for (i=0;i<numberExtraChangedBounds_;i++) {
1037    variables_[i]=rhs.variables_[i];
1038    newBounds_[i]=rhs.newBounds_[i];
1039  }
1040#endif
1041}
1042
1043// Assignment operator
1044CbcIntegerBranchingObject & 
1045CbcIntegerBranchingObject::operator=( const CbcIntegerBranchingObject& rhs)
1046{
1047  if (this != &rhs) {
1048    CbcBranchingObject::operator=(rhs);
1049    down_[0] = rhs.down_[0];
1050    down_[1] = rhs.down_[1];
1051    up_[0] = rhs.up_[0];
1052    up_[1] = rhs.up_[1];
1053#ifdef FUNNY_BRANCHING
1054    delete [] newBounds_;
1055    numberExtraChangedBounds_ = rhs.numberExtraChangedBounds_;
1056    int size = numberExtraChangedBounds_*(sizeof(double)+sizeof(int));
1057    char * temp = new char [size];
1058    newBounds_ = (double *) temp;
1059    variables_ = (int *) (newBounds_+numberExtraChangedBounds_);
1060   
1061    int i ;
1062    for (i=0;i<numberExtraChangedBounds_;i++) {
1063      variables_[i]=rhs.variables_[i];
1064      newBounds_[i]=rhs.newBounds_[i];
1065    }
1066#endif
1067  }
1068  return *this;
1069}
1070CbcBranchingObject * 
1071CbcIntegerBranchingObject::clone() const
1072{ 
1073  return (new CbcIntegerBranchingObject(*this));
1074}
1075
1076
1077// Destructor
1078CbcIntegerBranchingObject::~CbcIntegerBranchingObject ()
1079{
1080  // for debugging threads
1081  way_=-23456789;
1082#ifdef FUNNY_BRANCHING
1083  delete [] newBounds_;
1084#endif
1085}
1086
1087/*
1088  Perform a branch by adjusting the bounds of the specified variable. Note
1089  that each arm of the branch advances the object to the next arm by
1090  advancing the value of way_.
1091
1092  Providing new values for the variable's lower and upper bounds for each
1093  branching direction gives a little bit of additional flexibility and will
1094  be easily extensible to multi-way branching.
1095  Returns change in guessed objective on next branch
1096*/
1097double
1098CbcIntegerBranchingObject::branch()
1099{
1100  // for debugging threads
1101  if (way_<-1||way_>100000) {
1102    printf("way %d, left %d, iCol %d, variable %d\n",
1103           way_,numberBranchesLeft(),
1104           originalCbcObject_->columnNumber(),variable_);
1105    assert (way_!=-23456789);
1106  }
1107  decrementNumberBranchesLeft();
1108  if (down_[1]==-COIN_DBL_MAX)
1109    return 0.0;
1110  int iColumn = originalCbcObject_->columnNumber();
1111  assert (variable_==iColumn);
1112  double olb,oub ;
1113  olb = model_->solver()->getColLower()[iColumn] ;
1114  oub = model_->solver()->getColUpper()[iColumn] ;
1115#ifdef COIN_DEVELOP
1116  if (olb!=down_[0]||oub!=up_[1]) {
1117    if (way_>0)
1118      printf("branching up on var %d: [%g,%g] => [%g,%g] - other [%g,%g]\n",
1119             iColumn,olb,oub,up_[0],up_[1],down_[0],down_[1]) ; 
1120    else
1121      printf("branching down on var %d: [%g,%g] => [%g,%g] - other [%g,%g]\n",
1122             iColumn,olb,oub,down_[0],down_[1],up_[0],up_[1]) ; 
1123  }
1124#endif
1125  if (way_<0) {
1126#ifdef CBC_DEBUG
1127  { double olb,oub ;
1128    olb = model_->solver()->getColLower()[iColumn] ;
1129    oub = model_->solver()->getColUpper()[iColumn] ;
1130    printf("branching down on var %d: [%g,%g] => [%g,%g]\n",
1131           iColumn,olb,oub,down_[0],down_[1]) ; }
1132#endif
1133    model_->solver()->setColLower(iColumn,down_[0]);
1134    model_->solver()->setColUpper(iColumn,down_[1]);
1135    //#define CBC_PRINT2
1136#ifdef CBC_PRINT2
1137    printf("%d branching down has bounds %g %g",iColumn,down_[0],down_[1]);
1138#endif
1139#ifdef FUNNY_BRANCHING
1140    // branch - do extra bounds
1141    for (int i=0;i<numberExtraChangedBounds_;i++) {
1142      int variable = variables_[i];
1143      if ((variable&0x40000000)!=0) {
1144        // for going down
1145        int k = variable&0x3fffffff;
1146        assert (k!=iColumn);
1147        if ((variable&0x80000000)==0) {
1148          // lower bound changing
1149#ifdef CBC_PRINT2
1150          printf(" extra for %d changes lower from %g to %g",
1151                 k,model_->solver()->getColLower()[k],newBounds_[i]);
1152#endif
1153          model_->solver()->setColLower(k,newBounds_[i]);
1154        } else {
1155          // upper bound changing
1156#ifdef CBC_PRINT2
1157          printf(" extra for %d changes upper from %g to %g",
1158                 k,model_->solver()->getColUpper()[k],newBounds_[i]);
1159#endif
1160          model_->solver()->setColUpper(k,newBounds_[i]);
1161        }
1162      }
1163    }
1164#endif
1165#ifdef CBC_PRINT2
1166    printf("\n");
1167#endif
1168    way_=1;
1169  } else {
1170#ifdef CBC_DEBUG
1171  { double olb,oub ;
1172    olb = model_->solver()->getColLower()[iColumn] ;
1173    oub = model_->solver()->getColUpper()[iColumn] ;
1174    printf("branching up on var %d: [%g,%g] => [%g,%g]\n",
1175           iColumn,olb,oub,up_[0],up_[1]) ; }
1176#endif
1177    model_->solver()->setColLower(iColumn,up_[0]);
1178    model_->solver()->setColUpper(iColumn,up_[1]);
1179#ifdef CBC_PRINT2
1180    printf("%d branching up has bounds %g %g",iColumn,up_[0],up_[1]);
1181#endif
1182#ifdef FUNNY_BRANCHING
1183    // branch - do extra bounds
1184    for (int i=0;i<numberExtraChangedBounds_;i++) {
1185      int variable = variables_[i];
1186      if ((variable&0x40000000)==0) {
1187        // for going up
1188        int k = variable&0x3fffffff;
1189        assert (k!=iColumn);
1190        if ((variable&0x80000000)==0) {
1191          // lower bound changing
1192#ifdef CBC_PRINT2
1193          printf(" extra for %d changes lower from %g to %g",
1194                 k,model_->solver()->getColLower()[k],newBounds_[i]);
1195#endif
1196          model_->solver()->setColLower(k,newBounds_[i]);
1197        } else {
1198          // upper bound changing
1199#ifdef CBC_PRINT2
1200          printf(" extra for %d changes upper from %g to %g",
1201                 k,model_->solver()->getColUpper()[k],newBounds_[i]);
1202#endif
1203          model_->solver()->setColUpper(k,newBounds_[i]);
1204        }
1205      }
1206    }
1207#endif
1208#ifdef CBC_PRINT2
1209    printf("\n");
1210#endif
1211    way_=-1;      // Swap direction
1212  }
1213  double nlb = model_->solver()->getColLower()[iColumn];
1214  double nub = model_->solver()->getColUpper()[iColumn];
1215  if (nlb<olb) {
1216#ifndef NDEBUG
1217    printf("bad lb change for column %d from %g to %g\n",iColumn,olb,nlb);
1218#endif
1219    model_->solver()->setColLower(iColumn,CoinMin(olb,nub));
1220    nlb=olb;
1221  }
1222  if (nub>oub) {
1223#ifndef NDEBUG
1224    printf("bad ub change for column %d from %g to %g\n",iColumn,oub,nub);
1225#endif
1226    model_->solver()->setColUpper(iColumn,CoinMax(oub,nlb));
1227  }
1228#ifndef NDEBUG
1229  if (nlb<olb+1.0e-8&&nub>oub-1.0e-8&&false)
1230    printf("bad null change for column %d - bounds %g,%g\n",iColumn,olb,oub);
1231#endif
1232  return 0.0;
1233}
1234#ifdef FUNNY_BRANCHING
1235// Deactivate bounds for branching
1236void 
1237CbcIntegerBranchingObject::deactivate()
1238{
1239  down_[1]=-COIN_DBL_MAX;
1240}
1241int
1242CbcIntegerBranchingObject::applyExtraBounds(int iColumn, double lower, double upper, int way)
1243{
1244  // branch - do bounds
1245
1246  int i;
1247  int found=0;
1248  if (variable_==iColumn) {
1249    printf("odd applyExtra %d\n",iColumn);
1250    if (way<0) {
1251      down_[0]=CoinMax(lower,down_[0]);
1252      down_[1]=CoinMin(upper,down_[1]);
1253      assert (down_[0]<=down_[1]);
1254    } else {
1255      up_[0]=CoinMax(lower,up_[0]);
1256      up_[1]=CoinMin(upper,up_[1]);
1257      assert (up_[0]<=up_[1]);
1258    }
1259    return 0;
1260  }
1261  int check = (way<0) ? 0x40000000 : 0;
1262  double newLower=lower;
1263  double newUpper=upper;
1264  for (i=0;i<numberExtraChangedBounds_;i++) {
1265    int variable = variables_[i];
1266    if ((variable&0x40000000)==check) {
1267      int k = variable&0x3fffffff;
1268      if (k==iColumn) {
1269        if ((variable&0x80000000)==0) {
1270          // lower bound changing
1271          found |= 1;
1272          newBounds_[i] = CoinMax(lower,newBounds_[i]);
1273          newLower = newBounds_[i];
1274        } else {
1275          // upper bound changing
1276          found |= 2;
1277          newBounds_[i] = CoinMin(upper,newBounds_[i]);
1278          newUpper = newBounds_[i];
1279        }
1280      }
1281    }
1282  }
1283  int nAdd=0;
1284  if ((found&2)==0) {
1285    // need to add new upper
1286    nAdd++;
1287  }
1288  if ((found&1)==0) {
1289    // need to add new lower
1290    nAdd++;
1291  }
1292  if (nAdd) { 
1293    int size = (numberExtraChangedBounds_+nAdd)*(sizeof(double)+sizeof(int));
1294    char * temp = new char [size];
1295    double * newBounds = (double *) temp;
1296    int * variables = (int *) (newBounds+numberExtraChangedBounds_+nAdd);
1297
1298    int i ;
1299    for (i=0;i<numberExtraChangedBounds_;i++) {
1300      variables[i]=variables_[i];
1301      newBounds[i]=newBounds_[i];
1302    }
1303    delete [] newBounds_;
1304    newBounds_ = newBounds;
1305    variables_ = variables;
1306    if ((found&2)==0) {
1307      // need to add new upper
1308      int variable = iColumn | 0x80000000;
1309      variables_[numberExtraChangedBounds_]=variable;
1310      newBounds_[numberExtraChangedBounds_++]=newUpper;
1311    }
1312    if ((found&1)==0) {
1313      // need to add new lower
1314      int variable = iColumn;
1315      variables_[numberExtraChangedBounds_]=variable;
1316      newBounds_[numberExtraChangedBounds_++]=newLower;
1317    }
1318  }
1319 
1320  return (newUpper>=newLower) ? 0 : 1;
1321}
1322#endif
1323// Print what would happen 
1324void
1325CbcIntegerBranchingObject::print()
1326{
1327  int iColumn = originalCbcObject_->columnNumber();
1328  assert (variable_==iColumn);
1329  if (way_<0) {
1330  { double olb,oub ;
1331    olb = model_->solver()->getColLower()[iColumn] ;
1332    oub = model_->solver()->getColUpper()[iColumn] ;
1333    printf("CbcInteger would branch down on var %d (int var %d): [%g,%g] => [%g,%g]\n",
1334           iColumn,variable_,olb,oub,down_[0],down_[1]) ; }
1335  } else {
1336  { double olb,oub ;
1337    olb = model_->solver()->getColLower()[iColumn] ;
1338    oub = model_->solver()->getColUpper()[iColumn] ;
1339    printf("CbcInteger would branch up on var %d (int var %d): [%g,%g] => [%g,%g]\n",
1340           iColumn,variable_,olb,oub,up_[0],up_[1]) ; }
1341  }
1342}
1343
1344/** Compare the \c this with \c brObj. \c this and \c brObj must be os the
1345    same type and must have the same original object, but they may have
1346    different feasible regions.
1347    Return the appropriate CbcRangeCompare value (first argument being the
1348    sub/superset if that's the case). In case of overlap (and if \c
1349    replaceIfOverlap is true) replace the current branching object with one
1350    whose feasible region is the overlap.
1351   */
1352CbcRangeCompare
1353CbcIntegerBranchingObject::compareBranchingObject
1354(const CbcBranchingObject* brObj, const bool replaceIfOverlap)
1355{
1356  const CbcIntegerBranchingObject* br =
1357    dynamic_cast<const CbcIntegerBranchingObject*>(brObj);
1358  assert(br);
1359  double* thisBd = way_ < 0 ? down_ : up_;
1360  const double* otherBd = br->way_ < 0 ? br->down_ : br->up_;
1361  return CbcCompareRanges(thisBd, otherBd, replaceIfOverlap);
1362}
1363
1364//##############################################################################
1365
1366/** Default Constructor
1367
1368  Equivalent to an unspecified binary variable.
1369*/
1370CbcSimpleIntegerPseudoCost::CbcSimpleIntegerPseudoCost ()
1371  : CbcSimpleInteger(),
1372    downPseudoCost_(1.0e-5),
1373    upPseudoCost_(1.0e-5),
1374    upDownSeparator_(-1.0),
1375    method_(0)
1376{
1377}
1378
1379/** Useful constructor
1380
1381  Loads actual upper & lower bounds for the specified variable.
1382*/
1383CbcSimpleIntegerPseudoCost::CbcSimpleIntegerPseudoCost (CbcModel * model,
1384                                    int iColumn, double breakEven)
1385  : CbcSimpleInteger(model,iColumn,breakEven)
1386{
1387  const double * cost = model->getObjCoefficients();
1388  double costValue = CoinMax(1.0e-5,fabs(cost[iColumn]));
1389  // treat as if will cost what it says up
1390  upPseudoCost_=costValue;
1391  // and balance at breakeven
1392  downPseudoCost_=((1.0-breakEven_)*upPseudoCost_)/breakEven_;
1393  upDownSeparator_ = -1.0;
1394  method_=0;
1395}
1396
1397/** Useful constructor
1398
1399  Loads actual upper & lower bounds for the specified variable.
1400*/
1401CbcSimpleIntegerPseudoCost::CbcSimpleIntegerPseudoCost (CbcModel * model,
1402                                    int iColumn, double downPseudoCost,
1403                                                        double upPseudoCost)
1404  : CbcSimpleInteger(model,iColumn)
1405{
1406  downPseudoCost_ = CoinMax(1.0e-10,downPseudoCost);
1407  upPseudoCost_ = CoinMax(1.0e-10,upPseudoCost);
1408  breakEven_ = upPseudoCost_/(upPseudoCost_+downPseudoCost_);
1409  upDownSeparator_ = -1.0;
1410  method_=0;
1411}
1412// Useful constructor - passed and model index and pseudo costs
1413CbcSimpleIntegerPseudoCost::CbcSimpleIntegerPseudoCost (CbcModel * model, int dummy,int iColumn, 
1414                                                        double downPseudoCost, double upPseudoCost)
1415{
1416  *this=CbcSimpleIntegerPseudoCost(model,iColumn,downPseudoCost,upPseudoCost);
1417  columnNumber_=iColumn;
1418}
1419
1420// Copy constructor
1421CbcSimpleIntegerPseudoCost::CbcSimpleIntegerPseudoCost ( const CbcSimpleIntegerPseudoCost & rhs)
1422  :CbcSimpleInteger(rhs),
1423   downPseudoCost_(rhs.downPseudoCost_),
1424   upPseudoCost_(rhs.upPseudoCost_),
1425   upDownSeparator_(rhs.upDownSeparator_),
1426   method_(rhs.method_)
1427
1428{
1429}
1430
1431// Clone
1432CbcObject *
1433CbcSimpleIntegerPseudoCost::clone() const
1434{
1435  return new CbcSimpleIntegerPseudoCost(*this);
1436}
1437
1438// Assignment operator
1439CbcSimpleIntegerPseudoCost & 
1440CbcSimpleIntegerPseudoCost::operator=( const CbcSimpleIntegerPseudoCost& rhs)
1441{
1442  if (this!=&rhs) {
1443    CbcSimpleInteger::operator=(rhs);
1444    downPseudoCost_=rhs.downPseudoCost_;
1445    upPseudoCost_=rhs.upPseudoCost_;
1446    upDownSeparator_=rhs.upDownSeparator_;
1447    method_=rhs.method_;
1448  }
1449  return *this;
1450}
1451
1452// Destructor
1453CbcSimpleIntegerPseudoCost::~CbcSimpleIntegerPseudoCost ()
1454{
1455}
1456// Creates a branching object
1457CbcBranchingObject * 
1458CbcSimpleIntegerPseudoCost::createBranch(int way) 
1459{
1460  OsiSolverInterface * solver = model_->solver();
1461  const double * solution = model_->testSolution();
1462  const double * lower = solver->getColLower();
1463  const double * upper = solver->getColUpper();
1464  double value = solution[columnNumber_];
1465  value = CoinMax(value, lower[columnNumber_]);
1466  value = CoinMin(value, upper[columnNumber_]);
1467#ifndef NDEBUG
1468  double nearest = floor(value+0.5);
1469  double integerTolerance = 
1470    model_->getDblParam(CbcModel::CbcIntegerTolerance);
1471  assert (upper[columnNumber_]>lower[columnNumber_]);
1472#endif
1473  if (!model_->hotstartSolution()) {
1474    assert (fabs(value-nearest)>integerTolerance);
1475  } else {
1476    const double * hotstartSolution = model_->hotstartSolution();
1477    double targetValue = hotstartSolution[columnNumber_];
1478    if (way>0)
1479      value = targetValue-0.1;
1480    else
1481      value = targetValue+0.1;
1482  }
1483  CbcIntegerPseudoCostBranchingObject * newObject = 
1484    new CbcIntegerPseudoCostBranchingObject(model_,columnNumber_,way,
1485                                            value);
1486  double up =  upPseudoCost_*(ceil(value)-value);
1487  double down =  downPseudoCost_*(value-floor(value));
1488  double changeInGuessed=up-down;
1489  if (way>0)
1490    changeInGuessed = - changeInGuessed;
1491  changeInGuessed=CoinMax(0.0,changeInGuessed);
1492  //if (way>0)
1493  //changeInGuessed += 1.0e8; // bias to stay up
1494  newObject->setChangeInGuessed(changeInGuessed);
1495  newObject->setOriginalObject(this);
1496  return newObject;
1497}
1498// Infeasibility - large is 0.5
1499double 
1500CbcSimpleIntegerPseudoCost::infeasibility(int & preferredWay) const
1501{
1502  OsiSolverInterface * solver = model_->solver();
1503  const double * solution = model_->testSolution();
1504  const double * lower = solver->getColLower();
1505  const double * upper = solver->getColUpper();
1506  if (upper[columnNumber_]==lower[columnNumber_]) {
1507    // fixed
1508    preferredWay=1;
1509    return 0.0;
1510  }
1511  double value = solution[columnNumber_];
1512  value = CoinMax(value, lower[columnNumber_]);
1513  value = CoinMin(value, upper[columnNumber_]);
1514  /*printf("%d %g %g %g %g\n",columnNumber_,value,lower[columnNumber_],
1515    solution[columnNumber_],upper[columnNumber_]);*/
1516  double nearest = floor(value+0.5);
1517  double integerTolerance = 
1518    model_->getDblParam(CbcModel::CbcIntegerTolerance);
1519  double below = floor(value+integerTolerance);
1520  double above = below+1.0;
1521  if (above>upper[columnNumber_]) {
1522    above=below;
1523    below = above -1;
1524  }
1525  double downCost = CoinMax((value-below)*downPseudoCost_,0.0);
1526  double upCost = CoinMax((above-value)*upPseudoCost_,0.0);
1527  if (downCost>=upCost)
1528    preferredWay=1;
1529  else
1530    preferredWay=-1;
1531  // See if up down choice set
1532  if (upDownSeparator_>0.0) {
1533    preferredWay = (value-below>=upDownSeparator_) ? 1 : -1;
1534  }
1535  if (preferredWay_)
1536    preferredWay=preferredWay_;
1537  if (fabs(value-nearest)<=integerTolerance) {
1538    return 0.0;
1539  } else {
1540    // can't get at model so 1,2 don't make sense
1541    assert(method_<1||method_>2);
1542    if (!method_)
1543      return CoinMin(downCost,upCost);
1544    else
1545      return CoinMax(downCost,upCost);
1546  }
1547}
1548
1549// Return "up" estimate
1550double 
1551CbcSimpleIntegerPseudoCost::upEstimate() const
1552{
1553  OsiSolverInterface * solver = model_->solver();
1554  const double * solution = model_->testSolution();
1555  const double * lower = solver->getColLower();
1556  const double * upper = solver->getColUpper();
1557  double value = solution[columnNumber_];
1558  value = CoinMax(value, lower[columnNumber_]);
1559  value = CoinMin(value, upper[columnNumber_]);
1560  if (upper[columnNumber_]==lower[columnNumber_]) {
1561    // fixed
1562    return 0.0;
1563  }
1564  double integerTolerance = 
1565    model_->getDblParam(CbcModel::CbcIntegerTolerance);
1566  double below = floor(value+integerTolerance);
1567  double above = below+1.0;
1568  if (above>upper[columnNumber_]) {
1569    above=below;
1570    below = above -1;
1571  }
1572  double upCost = CoinMax((above-value)*upPseudoCost_,0.0);
1573  return upCost;
1574}
1575// Return "down" estimate
1576double 
1577CbcSimpleIntegerPseudoCost::downEstimate() const
1578{
1579  OsiSolverInterface * solver = model_->solver();
1580  const double * solution = model_->testSolution();
1581  const double * lower = solver->getColLower();
1582  const double * upper = solver->getColUpper();
1583  double value = solution[columnNumber_];
1584  value = CoinMax(value, lower[columnNumber_]);
1585  value = CoinMin(value, upper[columnNumber_]);
1586  if (upper[columnNumber_]==lower[columnNumber_]) {
1587    // fixed
1588    return 0.0;
1589  }
1590  double integerTolerance = 
1591    model_->getDblParam(CbcModel::CbcIntegerTolerance);
1592  double below = floor(value+integerTolerance);
1593  double above = below+1.0;
1594  if (above>upper[columnNumber_]) {
1595    above=below;
1596    below = above -1;
1597  }
1598  double downCost = CoinMax((value-below)*downPseudoCost_,0.0);
1599  return downCost;
1600}
1601
1602//##############################################################################
1603
1604// Default Constructor
1605CbcIntegerPseudoCostBranchingObject::CbcIntegerPseudoCostBranchingObject()
1606  :CbcIntegerBranchingObject()
1607{
1608  changeInGuessed_=1.0e-5;
1609}
1610
1611// Useful constructor
1612CbcIntegerPseudoCostBranchingObject::CbcIntegerPseudoCostBranchingObject (CbcModel * model, 
1613                                                      int variable, int way , double value)
1614  :CbcIntegerBranchingObject(model,variable,way,value)
1615{
1616}
1617// Useful constructor for fixing
1618CbcIntegerPseudoCostBranchingObject::CbcIntegerPseudoCostBranchingObject (CbcModel * model, 
1619                                                      int variable, int way,
1620                                                      double lowerValue, 
1621                                                      double upperValue)
1622  :CbcIntegerBranchingObject(model,variable,way,lowerValue)
1623{
1624  changeInGuessed_=1.0e100;
1625}
1626 
1627
1628// Copy constructor
1629CbcIntegerPseudoCostBranchingObject::CbcIntegerPseudoCostBranchingObject ( 
1630                                 const CbcIntegerPseudoCostBranchingObject & rhs)
1631  :CbcIntegerBranchingObject(rhs)
1632{
1633  changeInGuessed_ = rhs.changeInGuessed_;
1634}
1635
1636// Assignment operator
1637CbcIntegerPseudoCostBranchingObject & 
1638CbcIntegerPseudoCostBranchingObject::operator=( const CbcIntegerPseudoCostBranchingObject& rhs)
1639{
1640  if (this != &rhs) {
1641    CbcIntegerBranchingObject::operator=(rhs);
1642    changeInGuessed_ = rhs.changeInGuessed_;
1643  }
1644  return *this;
1645}
1646CbcBranchingObject * 
1647CbcIntegerPseudoCostBranchingObject::clone() const
1648{ 
1649  return (new CbcIntegerPseudoCostBranchingObject(*this));
1650}
1651
1652
1653// Destructor
1654CbcIntegerPseudoCostBranchingObject::~CbcIntegerPseudoCostBranchingObject ()
1655{
1656}
1657
1658/*
1659  Perform a branch by adjusting the bounds of the specified variable. Note
1660  that each arm of the branch advances the object to the next arm by
1661  advancing the value of way_.
1662
1663  Providing new values for the variable's lower and upper bounds for each
1664  branching direction gives a little bit of additional flexibility and will
1665  be easily extensible to multi-way branching.
1666  Returns change in guessed objective on next branch
1667*/
1668double
1669CbcIntegerPseudoCostBranchingObject::branch()
1670{
1671  CbcIntegerBranchingObject::branch();
1672  return changeInGuessed_;
1673}
1674
1675/** Compare the \c this with \c brObj. \c this and \c brObj must be os the
1676    same type and must have the same original object, but they may have
1677    different feasible regions.
1678    Return the appropriate CbcRangeCompare value (first argument being the
1679    sub/superset if that's the case). In case of overlap (and if \c
1680    replaceIfOverlap is true) replace the current branching object with one
1681    whose feasible region is the overlap.
1682*/
1683CbcRangeCompare
1684CbcIntegerPseudoCostBranchingObject::compareBranchingObject
1685(const CbcBranchingObject* brObj, const bool replaceIfOverlap)
1686{
1687  const CbcIntegerPseudoCostBranchingObject* br =
1688    dynamic_cast<const CbcIntegerPseudoCostBranchingObject*>(brObj);
1689  assert(br);
1690  double* thisBd = way_ < 0 ? down_ : up_;
1691  const double* otherBd = br->way_ < 0 ? br->down_ : br->up_;
1692  return CbcCompareRanges(thisBd, otherBd, replaceIfOverlap);
1693}
1694
1695
1696//##############################################################################
1697
1698// Default Constructor
1699CbcCliqueBranchingObject::CbcCliqueBranchingObject()
1700  :CbcBranchingObject()
1701{
1702  clique_ = NULL;
1703  downMask_[0]=0;
1704  downMask_[1]=0;
1705  upMask_[0]=0;
1706  upMask_[1]=0;
1707}
1708
1709// Useful constructor
1710CbcCliqueBranchingObject::CbcCliqueBranchingObject (CbcModel * model,
1711                                                    const CbcClique * clique,
1712                                                    int way ,
1713                                                    int numberOnDownSide, const int * down,
1714                                                    int numberOnUpSide, const int * up)
1715  :CbcBranchingObject(model,clique->id(),way,0.5)
1716{
1717  clique_ = clique;
1718  downMask_[0]=0;
1719  downMask_[1]=0;
1720  upMask_[0]=0;
1721  upMask_[1]=0;
1722  int i;
1723  for (i=0;i<numberOnDownSide;i++) {
1724    int sequence = down[i];
1725    int iWord = sequence>>5;
1726    int iBit = sequence - 32*iWord;
1727    unsigned int k = 1<<iBit;
1728    downMask_[iWord] |= k;
1729  }
1730  for (i=0;i<numberOnUpSide;i++) {
1731    int sequence = up[i];
1732    int iWord = sequence>>5;
1733    int iBit = sequence - 32*iWord;
1734    unsigned int k = 1<<iBit;
1735    upMask_[iWord] |= k;
1736  }
1737}
1738
1739// Copy constructor
1740CbcCliqueBranchingObject::CbcCliqueBranchingObject ( const CbcCliqueBranchingObject & rhs) :CbcBranchingObject(rhs)
1741{
1742  clique_=rhs.clique_;
1743  downMask_[0]=rhs.downMask_[0];
1744  downMask_[1]=rhs.downMask_[1];
1745  upMask_[0]=rhs.upMask_[0];
1746  upMask_[1]=rhs.upMask_[1];
1747}
1748
1749// Assignment operator
1750CbcCliqueBranchingObject & 
1751CbcCliqueBranchingObject::operator=( const CbcCliqueBranchingObject& rhs)
1752{
1753  if (this != &rhs) {
1754    CbcBranchingObject::operator=(rhs);
1755    clique_=rhs.clique_;
1756    downMask_[0]=rhs.downMask_[0];
1757    downMask_[1]=rhs.downMask_[1];
1758    upMask_[0]=rhs.upMask_[0];
1759    upMask_[1]=rhs.upMask_[1];
1760  }
1761  return *this;
1762}
1763CbcBranchingObject * 
1764CbcCliqueBranchingObject::clone() const
1765{ 
1766  return (new CbcCliqueBranchingObject(*this));
1767}
1768
1769
1770// Destructor
1771CbcCliqueBranchingObject::~CbcCliqueBranchingObject ()
1772{
1773}
1774double
1775CbcCliqueBranchingObject::branch()
1776{
1777  decrementNumberBranchesLeft();
1778  int iWord;
1779  int numberMembers = clique_->numberMembers();
1780  const int * which = clique_->members();
1781  const int * integerVariables = model_->integerVariable();
1782  int numberWords=(numberMembers+31)>>5;
1783  // *** for way - up means fix all those in down section
1784  if (way_<0) {
1785#ifdef FULL_PRINT
1786    printf("Down Fix ");
1787#endif
1788    for (iWord=0;iWord<numberWords;iWord++) {
1789      int i;
1790      for (i=0;i<32;i++) {
1791        unsigned int k = 1<<i;
1792        if ((upMask_[iWord]&k)!=0) {
1793          int iColumn = which[i+32*iWord];
1794#ifdef FULL_PRINT
1795          printf("%d ",i+32*iWord);
1796#endif
1797          // fix weak way
1798          if (clique_->type(i+32*iWord))
1799            model_->solver()->setColUpper(integerVariables[iColumn],0.0);
1800          else
1801            model_->solver()->setColLower(integerVariables[iColumn],1.0);
1802        }
1803      }
1804    }
1805    way_=1;       // Swap direction
1806  } else {
1807#ifdef FULL_PRINT
1808    printf("Up Fix ");
1809#endif
1810    for (iWord=0;iWord<numberWords;iWord++) {
1811      int i;
1812      for (i=0;i<32;i++) {
1813        unsigned int k = 1<<i;
1814        if ((downMask_[iWord]&k)!=0) {
1815          int iColumn = which[i+32*iWord];
1816#ifdef FULL_PRINT
1817          printf("%d ",i+32*iWord);
1818#endif
1819          // fix weak way
1820          if (clique_->type(i+32*iWord))
1821            model_->solver()->setColUpper(integerVariables[iColumn],0.0);
1822          else
1823            model_->solver()->setColLower(integerVariables[iColumn],1.0);
1824        }
1825      }
1826    }
1827    way_=-1;      // Swap direction
1828  }
1829#ifdef FULL_PRINT
1830  printf("\n");
1831#endif
1832  return 0.0;
1833}
1834// Print what would happen 
1835void
1836CbcCliqueBranchingObject::print()
1837{
1838  int iWord;
1839  int numberMembers = clique_->numberMembers();
1840  const int * which = clique_->members();
1841  const int * integerVariables = model_->integerVariable();
1842  int numberWords=(numberMembers+31)>>5;
1843  // *** for way - up means fix all those in down section
1844  if (way_<0) {
1845    printf("Clique - Down Fix ");
1846    for (iWord=0;iWord<numberWords;iWord++) {
1847      int i;
1848      for (i=0;i<32;i++) {
1849        unsigned int k = 1<<i;
1850        if ((upMask_[iWord]&k)!=0) {
1851          int iColumn = which[i+32*iWord];
1852          printf("%d ",integerVariables[iColumn]);
1853        }
1854      }
1855    }
1856  } else {
1857    printf("Clique - Up Fix ");
1858    for (iWord=0;iWord<numberWords;iWord++) {
1859      int i;
1860      for (i=0;i<32;i++) {
1861        unsigned int k = 1<<i;
1862        if ((downMask_[iWord]&k)!=0) {
1863          int iColumn = which[i+32*iWord];
1864          printf("%d ",integerVariables[iColumn]);
1865        }
1866      }
1867    }
1868  }
1869  printf("\n");
1870}
1871
1872static inline int
1873CbcCompareCliques(const CbcClique* cl0, const CbcClique* cl1)
1874{
1875  if (cl0->cliqueType() < cl1->cliqueType()) {
1876    return -1;
1877  }
1878  if (cl0->cliqueType() > cl1->cliqueType()) {
1879    return 1;
1880  }
1881  if (cl0->numberMembers() != cl1->numberMembers()) {
1882    return cl0->numberMembers() - cl1->numberMembers();
1883  }
1884  if (cl0->numberNonSOSMembers() != cl1->numberNonSOSMembers()) {
1885    return cl0->numberNonSOSMembers() - cl1->numberNonSOSMembers();
1886  }
1887  return memcmp(cl0->members(), cl1->members(),
1888                cl0->numberMembers() * sizeof(int));
1889}
1890
1891/** Compare the original object of \c this with the original object of \c
1892    brObj. Assumes that there is an ordering of the original objects.
1893    This method should be invoked only if \c this and brObj are of the same
1894    type.
1895    Return negative/0/positive depending on whether \c this is
1896    smaller/same/larger than the argument.
1897*/
1898int
1899CbcCliqueBranchingObject::compareOriginalObject
1900(const CbcBranchingObject* brObj) const
1901{
1902  const CbcCliqueBranchingObject* br =
1903    dynamic_cast<const CbcCliqueBranchingObject*>(brObj);
1904  assert(br);
1905  return CbcCompareCliques(clique_, br->clique_);
1906}
1907
1908/** Compare the \c this with \c brObj. \c this and \c brObj must be os the
1909    same type and must have the same original object, but they may have
1910    different feasible regions.
1911    Return the appropriate CbcRangeCompare value (first argument being the
1912    sub/superset if that's the case). In case of overlap (and if \c
1913    replaceIfOverlap is true) replace the current branching object with one
1914    whose feasible region is the overlap.
1915*/
1916CbcRangeCompare
1917CbcCliqueBranchingObject::compareBranchingObject
1918(const CbcBranchingObject* brObj, const bool replaceIfOverlap)
1919{
1920  const CbcCliqueBranchingObject* br =
1921    dynamic_cast<const CbcCliqueBranchingObject*>(brObj);
1922  assert(br);
1923  unsigned int* thisMask = way_ < 0 ? upMask_ : downMask_;
1924  const unsigned int* otherMask = br->way_ < 0 ? br->upMask_ : br->downMask_;
1925  const CoinUInt64 cl0 =
1926    (static_cast<CoinUInt64>(thisMask[0]) << 32) | thisMask[1];
1927  const CoinUInt64 cl1 =
1928    (static_cast<CoinUInt64>(otherMask[0]) << 32) | otherMask[1];
1929  if (cl0 == cl1) {
1930    return CbcRangeSame;
1931  }
1932  const CoinUInt64 cl_intersection = (cl0 & cl1);
1933  if (cl_intersection == cl0) {
1934    return CbcRangeSuperset;
1935  }
1936  if (cl_intersection == cl1) {
1937    return CbcRangeSubset;
1938  }
1939  const CoinUInt64 cl_xor = (cl0 ^ cl1);
1940  if (cl_intersection == 0 && cl_xor == 0) {
1941    return CbcRangeDisjoint;
1942  }
1943  const CoinUInt64 cl_union = (cl0 | cl1);
1944  thisMask[0] = static_cast<unsigned int>(cl_union >> 32);
1945  thisMask[1] = static_cast<unsigned int>(cl_union & 0xffffffff);
1946  return CbcRangeOverlap;
1947}
1948
1949//##############################################################################
1950
1951// Default Constructor
1952CbcLongCliqueBranchingObject::CbcLongCliqueBranchingObject()
1953  :CbcBranchingObject()
1954{
1955  clique_=NULL;
1956  downMask_=NULL;
1957  upMask_=NULL;
1958}
1959
1960// Useful constructor
1961CbcLongCliqueBranchingObject::CbcLongCliqueBranchingObject (CbcModel * model,
1962                                                            const CbcClique * clique, 
1963                                                            int way ,
1964                                                    int numberOnDownSide, const int * down,
1965                                                    int numberOnUpSide, const int * up)
1966  :CbcBranchingObject(model,clique->id(),way,0.5)
1967{
1968  clique_ = clique;
1969  int numberMembers = clique_->numberMembers();
1970  int numberWords=(numberMembers+31)>>5;
1971  downMask_ = new unsigned int [numberWords];
1972  upMask_ = new unsigned int [numberWords];
1973  memset(downMask_,0,numberWords*sizeof(unsigned int));
1974  memset(upMask_,0,numberWords*sizeof(unsigned int));
1975  int i;
1976  for (i=0;i<numberOnDownSide;i++) {
1977    int sequence = down[i];
1978    int iWord = sequence>>5;
1979    int iBit = sequence - 32*iWord;
1980    unsigned int k = 1<<iBit;
1981    downMask_[iWord] |= k;
1982  }
1983  for (i=0;i<numberOnUpSide;i++) {
1984    int sequence = up[i];
1985    int iWord = sequence>>5;
1986    int iBit = sequence - 32*iWord;
1987    unsigned int k = 1<<iBit;
1988    upMask_[iWord] |= k;
1989  }
1990}
1991
1992// Copy constructor
1993CbcLongCliqueBranchingObject::CbcLongCliqueBranchingObject ( const CbcLongCliqueBranchingObject & rhs) :CbcBranchingObject(rhs)
1994{
1995  clique_=rhs.clique_;
1996  if (rhs.downMask_) {
1997    int numberMembers = clique_->numberMembers();
1998    int numberWords=(numberMembers+31)>>5;
1999    downMask_ = new unsigned int [numberWords];
2000    memcpy(downMask_,rhs.downMask_,numberWords*sizeof(unsigned int));
2001    upMask_ = new unsigned int [numberWords];
2002    memcpy(upMask_,rhs.upMask_,numberWords*sizeof(unsigned int));
2003  } else {
2004    downMask_=NULL;
2005    upMask_=NULL;
2006  }   
2007}
2008
2009// Assignment operator
2010CbcLongCliqueBranchingObject & 
2011CbcLongCliqueBranchingObject::operator=( const CbcLongCliqueBranchingObject& rhs)
2012{
2013  if (this != &rhs) {
2014    CbcBranchingObject::operator=(rhs);
2015    clique_=rhs.clique_;
2016    delete [] downMask_;
2017    delete [] upMask_;
2018    if (rhs.downMask_) {
2019      int numberMembers = clique_->numberMembers();
2020      int numberWords=(numberMembers+31)>>5;
2021      downMask_ = new unsigned int [numberWords];
2022      memcpy(downMask_,rhs.downMask_,numberWords*sizeof(unsigned int));
2023      upMask_ = new unsigned int [numberWords];
2024      memcpy(upMask_,rhs.upMask_,numberWords*sizeof(unsigned int));
2025    } else {
2026      downMask_=NULL;
2027      upMask_=NULL;
2028    }   
2029  }
2030  return *this;
2031}
2032CbcBranchingObject * 
2033CbcLongCliqueBranchingObject::clone() const
2034{ 
2035  return (new CbcLongCliqueBranchingObject(*this));
2036}
2037
2038
2039// Destructor
2040CbcLongCliqueBranchingObject::~CbcLongCliqueBranchingObject ()
2041{
2042  delete [] downMask_;
2043  delete [] upMask_;
2044}
2045double
2046CbcLongCliqueBranchingObject::branch()
2047{
2048  decrementNumberBranchesLeft();
2049  int iWord;
2050  int numberMembers = clique_->numberMembers();
2051  const int * which = clique_->members();
2052  const int * integerVariables = model_->integerVariable();
2053  int numberWords=(numberMembers+31)>>5;
2054  // *** for way - up means fix all those in down section
2055  if (way_<0) {
2056#ifdef FULL_PRINT
2057    printf("Down Fix ");
2058#endif
2059    for (iWord=0;iWord<numberWords;iWord++) {
2060      int i;
2061      for (i=0;i<32;i++) {
2062        unsigned int k = 1<<i;
2063        if ((upMask_[iWord]&k)!=0) {
2064          int iColumn = which[i+32*iWord];
2065#ifdef FULL_PRINT
2066          printf("%d ",i+32*iWord);
2067#endif
2068          // fix weak way
2069          if (clique_->type(i+32*iWord))
2070            model_->solver()->setColUpper(integerVariables[iColumn],0.0);
2071          else
2072            model_->solver()->setColLower(integerVariables[iColumn],1.0);
2073        }
2074      }
2075    }
2076    way_=1;       // Swap direction
2077  } else {
2078#ifdef FULL_PRINT
2079    printf("Up Fix ");
2080#endif
2081    for (iWord=0;iWord<numberWords;iWord++) {
2082      int i;
2083      for (i=0;i<32;i++) {
2084        unsigned int k = 1<<i;
2085        if ((downMask_[iWord]&k)!=0) {
2086          int iColumn = which[i+32*iWord];
2087#ifdef FULL_PRINT
2088          printf("%d ",i+32*iWord);
2089#endif
2090          // fix weak way
2091          if (clique_->type(i+32*iWord))
2092            model_->solver()->setColUpper(integerVariables[iColumn],0.0);
2093          else
2094            model_->solver()->setColLower(integerVariables[iColumn],1.0);
2095        }
2096      }
2097    }
2098    way_=-1;      // Swap direction
2099  }
2100#ifdef FULL_PRINT
2101  printf("\n");
2102#endif
2103  return 0.0;
2104}
2105void
2106CbcLongCliqueBranchingObject::print()
2107{
2108  int iWord;
2109  int numberMembers = clique_->numberMembers();
2110  const int * which = clique_->members();
2111  const int * integerVariables = model_->integerVariable();
2112  int numberWords=(numberMembers+31)>>5;
2113  // *** for way - up means fix all those in down section
2114  if (way_<0) {
2115    printf("Clique - Down Fix ");
2116    for (iWord=0;iWord<numberWords;iWord++) {
2117      int i;
2118      for (i=0;i<32;i++) {
2119        unsigned int k = 1<<i;
2120        if ((upMask_[iWord]&k)!=0) {
2121          int iColumn = which[i+32*iWord];
2122          printf("%d ",integerVariables[iColumn]);
2123        }
2124      }
2125    }
2126  } else {
2127    printf("Clique - Up Fix ");
2128    for (iWord=0;iWord<numberWords;iWord++) {
2129      int i;
2130      for (i=0;i<32;i++) {
2131        unsigned int k = 1<<i;
2132        if ((downMask_[iWord]&k)!=0) {
2133          int iColumn = which[i+32*iWord];
2134          printf("%d ",integerVariables[iColumn]);
2135        }
2136      }
2137    }
2138  }
2139  printf("\n");
2140}
2141
2142/** Compare the original object of \c this with the original object of \c
2143    brObj. Assumes that there is an ordering of the original objects.
2144    This method should be invoked only if \c this and brObj are of the same
2145    type.
2146    Return negative/0/positive depending on whether \c this is
2147    smaller/same/larger than the argument.
2148*/
2149int
2150CbcLongCliqueBranchingObject::compareOriginalObject
2151(const CbcBranchingObject* brObj) const
2152{
2153  const CbcLongCliqueBranchingObject* br =
2154    dynamic_cast<const CbcLongCliqueBranchingObject*>(brObj);
2155  assert(br);
2156  return CbcCompareCliques(clique_, br->clique_);
2157}
2158
2159/** Compare the \c this with \c brObj. \c this and \c brObj must be os the
2160    same type and must have the same original object, but they may have
2161    different feasible regions.
2162    Return the appropriate CbcRangeCompare value (first argument being the
2163    sub/superset if that's the case). In case of overlap (and if \c
2164    replaceIfOverlap is true) replace the current branching object with one
2165    whose feasible region is the overlap.
2166*/
2167CbcRangeCompare
2168CbcLongCliqueBranchingObject::compareBranchingObject
2169(const CbcBranchingObject* brObj, const bool replaceIfOverlap)
2170{
2171  const CbcLongCliqueBranchingObject* br =
2172    dynamic_cast<const CbcLongCliqueBranchingObject*>(brObj);
2173  assert(br);
2174  const int numberMembers = clique_->numberMembers();
2175  const int numberWords=(numberMembers+31)>>5;
2176  unsigned int* thisMask = way_ < 0 ? upMask_ : downMask_;
2177  const unsigned int* otherMask = br->way_ < 0 ? br->upMask_ : br->downMask_;
2178
2179  if (memcmp(thisMask, otherMask, numberWords * sizeof(unsigned int)) == 0) {
2180    return CbcRangeSame;
2181  }
2182  bool canBeSuperset = true;
2183  bool canBeSubset = true;
2184  int i;
2185  for (i = numberWords-1; i >= 0 && (canBeSuperset || canBeSubset); --i) {
2186    const unsigned int both = (thisMask[i] & otherMask[i]);
2187    canBeSuperset &= (both == thisMask[i]);
2188    canBeSubset &= (both == otherMask[i]);
2189  }
2190  if (canBeSuperset) {
2191    return CbcRangeSuperset;
2192  }
2193  if (canBeSubset) {
2194    return CbcRangeSubset;
2195  }
2196
2197  for (i = numberWords-1; i >= 0; --i) {
2198    if ((thisMask[i] ^ otherMask[i]) != 0) {
2199      break;
2200    }
2201  }
2202  if (i == -1) { // complement
2203    return CbcRangeDisjoint;
2204  }
2205  // must be overlap
2206  for (i = numberWords-1; i >= 0; --i) {
2207    thisMask[i] |= otherMask[i];
2208  }
2209  return CbcRangeOverlap;
2210}
2211
2212//##############################################################################
2213
2214// Default Constructor
2215CbcSOSBranchingObject::CbcSOSBranchingObject()
2216  :CbcBranchingObject(),
2217   firstNonzero_(-1),
2218   lastNonzero_(-1)
2219{
2220  set_ = NULL;
2221  separator_=0.0;
2222}
2223
2224// Useful constructor
2225CbcSOSBranchingObject::CbcSOSBranchingObject (CbcModel * model,
2226                                              const CbcSOS * set,
2227                                              int way ,
2228                                              double separator)
2229  :CbcBranchingObject(model,set->id(),way,0.5)
2230{
2231  set_ = set;
2232  separator_ = separator;
2233  computeNonzeroRange();
2234}
2235
2236// Copy constructor
2237CbcSOSBranchingObject::CbcSOSBranchingObject (const CbcSOSBranchingObject & rhs)
2238 :CbcBranchingObject(rhs),
2239  firstNonzero_(rhs.firstNonzero_),
2240  lastNonzero_(rhs.lastNonzero_)
2241{
2242  set_=rhs.set_;
2243  separator_ = rhs.separator_;
2244}
2245
2246// Assignment operator
2247CbcSOSBranchingObject & 
2248CbcSOSBranchingObject::operator=( const CbcSOSBranchingObject& rhs)
2249{
2250  if (this != &rhs) {
2251    CbcBranchingObject::operator=(rhs);
2252    set_=rhs.set_;
2253    separator_ = rhs.separator_;
2254    firstNonzero_ = rhs.firstNonzero_;
2255    lastNonzero_ = rhs.lastNonzero_;
2256  }
2257  return *this;
2258}
2259CbcBranchingObject * 
2260CbcSOSBranchingObject::clone() const
2261{ 
2262  return (new CbcSOSBranchingObject(*this));
2263}
2264
2265
2266// Destructor
2267CbcSOSBranchingObject::~CbcSOSBranchingObject ()
2268{
2269}
2270
2271void
2272CbcSOSBranchingObject::computeNonzeroRange()
2273{
2274  const int numberMembers = set_->numberMembers();
2275  const double * weights = set_->weights();
2276  int i = 0;
2277  if (way_ < 0) {
2278    for ( i=0;i<numberMembers;i++) {
2279      if (weights[i] > separator_)
2280        break;
2281    }
2282    assert (i<numberMembers);
2283    firstNonzero_ = 0;
2284    lastNonzero_ = i;
2285  } else {
2286    for ( i=0;i<numberMembers;i++) {
2287      if (weights[i] >= separator_)
2288        break;
2289    }
2290    assert (i<numberMembers);
2291    firstNonzero_ = i;
2292    lastNonzero_ = numberMembers;
2293  }
2294}
2295
2296double
2297CbcSOSBranchingObject::branch()
2298{
2299  decrementNumberBranchesLeft();
2300  int numberMembers = set_->numberMembers();
2301  const int * which = set_->members();
2302  const double * weights = set_->weights();
2303  OsiSolverInterface * solver = model_->solver();
2304  //const double * lower = solver->getColLower();
2305  //const double * upper = solver->getColUpper();
2306  // *** for way - up means fix all those in down section
2307  if (way_<0) {
2308    int i;
2309    for ( i=0;i<numberMembers;i++) {
2310      if (weights[i] > separator_)
2311        break;
2312    }
2313    assert (i<numberMembers);
2314    for (;i<numberMembers;i++) 
2315      solver->setColUpper(which[i],0.0);
2316    way_=1;       // Swap direction
2317  } else {
2318    int i;
2319    for ( i=0;i<numberMembers;i++) {
2320      if (weights[i] >= separator_)
2321        break;
2322      else
2323        solver->setColUpper(which[i],0.0);
2324    }
2325    assert (i<numberMembers);
2326    way_=-1;      // Swap direction
2327  }
2328  computeNonzeroRange();
2329  return 0.0;
2330}
2331// Print what would happen 
2332void
2333CbcSOSBranchingObject::print()
2334{
2335  int numberMembers = set_->numberMembers();
2336  const int * which = set_->members();
2337  const double * weights = set_->weights();
2338  OsiSolverInterface * solver = model_->solver();
2339  //const double * lower = solver->getColLower();
2340  const double * upper = solver->getColUpper();
2341  int first=numberMembers;
2342  int last=-1;
2343  int numberFixed=0;
2344  int numberOther=0;
2345  int i;
2346  for ( i=0;i<numberMembers;i++) {
2347    double bound = upper[which[i]];
2348    if (bound) {
2349      first = CoinMin(first,i);
2350      last = CoinMax(last,i);
2351    }
2352  }
2353  // *** for way - up means fix all those in down section
2354  if (way_<0) {
2355    printf("SOS Down");
2356    for ( i=0;i<numberMembers;i++) {
2357      double bound = upper[which[i]];
2358      if (weights[i] > separator_)
2359        break;
2360      else if (bound)
2361        numberOther++;
2362    }
2363    assert (i<numberMembers);
2364    for (;i<numberMembers;i++) {
2365      double bound = upper[which[i]];
2366      if (bound)
2367        numberFixed++;
2368    }
2369  } else {
2370    printf("SOS Up");
2371    for ( i=0;i<numberMembers;i++) {
2372      double bound = upper[which[i]];
2373      if (weights[i] >= separator_)
2374        break;
2375      else if (bound)
2376        numberFixed++;
2377    }
2378    assert (i<numberMembers);
2379    for (;i<numberMembers;i++) {
2380      double bound = upper[which[i]];
2381      if (bound)
2382        numberOther++;
2383    }
2384  }
2385  printf(" - at %g, free range %d (%g) => %d (%g), %d would be fixed, %d other way\n",
2386         separator_,which[first],weights[first],which[last],weights[last],numberFixed,numberOther);
2387}
2388 
2389/** Compare the original object of \c this with the original object of \c
2390    brObj. Assumes that there is an ordering of the original objects.
2391    This method should be invoked only if \c this and brObj are of the same
2392    type.
2393    Return negative/0/positive depending on whether \c this is
2394    smaller/same/larger than the argument.
2395*/
2396int
2397CbcSOSBranchingObject::compareOriginalObject
2398(const CbcBranchingObject* brObj) const
2399{
2400  const CbcSOSBranchingObject* br =
2401    dynamic_cast<const CbcSOSBranchingObject*>(brObj);
2402  assert(br);
2403  const CbcSOS* s0 = set_;
2404  const CbcSOS* s1 = br->set_;
2405  if (s0->sosType() != s1->sosType()) {
2406    return s0->sosType() - s1->sosType();
2407  }
2408  if (s0->numberMembers() != s1->numberMembers()) {
2409    return s0->numberMembers() - s1->numberMembers();
2410  }
2411  const int memberCmp = memcmp(s0->members(), s1->members(),
2412                               s0->numberMembers() * sizeof(int));
2413  if (memberCmp != 0) {
2414    return memberCmp;
2415  }
2416  return memcmp(s0->weights(), s1->weights(),
2417                s0->numberMembers() * sizeof(double));
2418}
2419
2420/** Compare the \c this with \c brObj. \c this and \c brObj must be os the
2421    same type and must have the same original object, but they may have
2422    different feasible regions.
2423    Return the appropriate CbcRangeCompare value (first argument being the
2424    sub/superset if that's the case). In case of overlap (and if \c
2425    replaceIfOverlap is true) replace the current branching object with one
2426    whose feasible region is the overlap.
2427*/
2428CbcRangeCompare
2429CbcSOSBranchingObject::compareBranchingObject
2430(const CbcBranchingObject* brObj, const bool replaceIfOverlap)
2431{
2432  const CbcSOSBranchingObject* br =
2433    dynamic_cast<const CbcSOSBranchingObject*>(brObj);
2434  assert(br);
2435  if (firstNonzero_ < br->firstNonzero_) {
2436    if (lastNonzero_ >= br->lastNonzero_) {
2437      return CbcRangeSuperset;
2438    } else if (lastNonzero_ <= br->firstNonzero_) {
2439      return CbcRangeDisjoint;
2440    } else {
2441      // overlap
2442      if (replaceIfOverlap) {
2443        firstNonzero_ = br->firstNonzero_;
2444      }
2445      return CbcRangeOverlap;
2446    }
2447  } else if (firstNonzero_ > br->firstNonzero_) {
2448    if (lastNonzero_ <= br->lastNonzero_) {
2449      return CbcRangeSubset;
2450    } else if (firstNonzero_ >= br->lastNonzero_) {
2451      return CbcRangeDisjoint;
2452    } else {
2453      // overlap
2454      if (replaceIfOverlap) {
2455        lastNonzero_ = br->lastNonzero_;
2456      }
2457      return CbcRangeOverlap;
2458    }
2459  } else {
2460    if (lastNonzero_ == br->lastNonzero_) {
2461      return CbcRangeSame;
2462    }
2463    return lastNonzero_ < br->lastNonzero_ ? CbcRangeSubset : CbcRangeSuperset;
2464  }
2465  return CbcRangeSame; // fake return
2466}
2467
2468//##############################################################################
2469
2470// Default Constructor
2471CbcBranchDefaultDecision::CbcBranchDefaultDecision()
2472  :CbcBranchDecision()
2473{
2474  bestCriterion_ = 0.0;
2475  bestChangeUp_ = 0.0;
2476  bestNumberUp_ = 0;
2477  bestChangeDown_ = 0.0;
2478  bestObject_ = NULL;
2479  bestNumberDown_ = 0;
2480}
2481
2482// Copy constructor
2483CbcBranchDefaultDecision::CbcBranchDefaultDecision (
2484                                    const CbcBranchDefaultDecision & rhs)
2485  :CbcBranchDecision(rhs)
2486{
2487  bestCriterion_ = rhs.bestCriterion_;
2488  bestChangeUp_ = rhs.bestChangeUp_;
2489  bestNumberUp_ = rhs.bestNumberUp_;
2490  bestChangeDown_ = rhs.bestChangeDown_;
2491  bestNumberDown_ = rhs.bestNumberDown_;
2492  bestObject_ = rhs.bestObject_;
2493  model_ = rhs.model_;
2494}
2495
2496CbcBranchDefaultDecision::~CbcBranchDefaultDecision()
2497{
2498}
2499
2500// Clone
2501CbcBranchDecision * 
2502CbcBranchDefaultDecision::clone() const
2503{
2504  return new CbcBranchDefaultDecision(*this);
2505}
2506
2507// Initialize i.e. before start of choosing at a node
2508void 
2509CbcBranchDefaultDecision::initialize(CbcModel * model)
2510{
2511  bestCriterion_ = 0.0;
2512  bestChangeUp_ = 0.0;
2513  bestNumberUp_ = 0;
2514  bestChangeDown_ = 0.0;
2515  bestNumberDown_ = 0;
2516  bestObject_ = NULL;
2517  model_ = model;
2518}
2519
2520
2521/*
2522  Simple default decision algorithm. Compare based on infeasibility (numInfUp,
2523  numInfDn) until a solution is found by search, then switch to change in
2524  objective (changeUp, changeDn). Note that bestSoFar is remembered in
2525  bestObject_, so the parameter bestSoFar is unused.
2526*/
2527
2528int
2529CbcBranchDefaultDecision::betterBranch(CbcBranchingObject * thisOne,
2530                            CbcBranchingObject * bestSoFar,
2531                            double changeUp, int numInfUp,
2532                            double changeDn, int numInfDn)
2533{
2534  bool beforeSolution = cbcModel()->getSolutionCount()==
2535    cbcModel()->getNumberHeuristicSolutions();;
2536  int betterWay=0;
2537  if (beforeSolution) {
2538    if (!bestObject_) {
2539      bestNumberUp_=COIN_INT_MAX;
2540      bestNumberDown_=COIN_INT_MAX;
2541    }
2542    // before solution - choose smallest number
2543    // could add in depth as well
2544    int bestNumber = CoinMin(bestNumberUp_,bestNumberDown_);
2545    if (numInfUp<numInfDn) {
2546      if (numInfUp<bestNumber) {
2547        betterWay = 1;
2548      } else if (numInfUp==bestNumber) {
2549        if (changeUp<bestCriterion_)
2550          betterWay=1;
2551      }
2552    } else if (numInfUp>numInfDn) {
2553      if (numInfDn<bestNumber) {
2554        betterWay = -1;
2555      } else if (numInfDn==bestNumber) {
2556        if (changeDn<bestCriterion_)
2557          betterWay=-1;
2558      }
2559    } else {
2560      // up and down have same number
2561      bool better=false;
2562      if (numInfUp<bestNumber) {
2563        better=true;
2564      } else if (numInfUp==bestNumber) {
2565        if (min(changeUp,changeDn)<bestCriterion_)
2566          better=true;;
2567      }
2568      if (better) {
2569        // see which way
2570        if (changeUp<=changeDn)
2571          betterWay=1;
2572        else
2573          betterWay=-1;
2574      }
2575    }
2576  } else {
2577    if (!bestObject_) {
2578      bestCriterion_=-1.0;
2579    }
2580    // got a solution
2581    if (changeUp<=changeDn) {
2582      if (changeUp>bestCriterion_)
2583        betterWay=1;
2584    } else {
2585      if (changeDn>bestCriterion_)
2586        betterWay=-1;
2587    }
2588  }
2589  if (betterWay) {
2590    bestCriterion_ = CoinMin(changeUp,changeDn);
2591    bestChangeUp_ = changeUp;
2592    bestNumberUp_ = numInfUp;
2593    bestChangeDown_ = changeDn;
2594    bestNumberDown_ = numInfDn;
2595    bestObject_=thisOne;
2596    // See if user is overriding way
2597    if (thisOne->object()&&thisOne->object()->preferredWay())
2598      betterWay = thisOne->object()->preferredWay();
2599  }
2600  return betterWay;
2601}
2602/* Sets or gets best criterion so far */
2603void 
2604CbcBranchDefaultDecision::setBestCriterion(double value)
2605{ 
2606  bestCriterion_ = value;
2607}
2608double 
2609CbcBranchDefaultDecision::getBestCriterion() const
2610{ 
2611  return bestCriterion_;
2612}
2613
2614/* Compare N branching objects. Return index of best
2615   and sets way of branching in chosen object.
2616   
2617   This routine is used only after strong branching.
2618*/
2619
2620int
2621CbcBranchDefaultDecision::bestBranch (CbcBranchingObject ** objects, int numberObjects,
2622                                   int numberUnsatisfied,
2623                                   double * changeUp, int * numberInfeasibilitiesUp,
2624                                   double * changeDown, int * numberInfeasibilitiesDown,
2625                                   double objectiveValue) 
2626{
2627
2628  int bestWay=0;
2629  int whichObject = -1;
2630  if (numberObjects) {
2631    CbcModel * model = cbcModel();
2632    // at continuous
2633    //double continuousObjective = model->getContinuousObjective();
2634    //int continuousInfeasibilities = model->getContinuousInfeasibilities();
2635   
2636    // average cost to get rid of infeasibility
2637    //double averageCostPerInfeasibility =
2638    //(objectiveValue-continuousObjective)/
2639    //(double) (abs(continuousInfeasibilities-numberUnsatisfied)+1);
2640    /* beforeSolution is :
2641       0 - before any solution
2642       n - n heuristic solutions but no branched one
2643       -1 - branched solution found
2644    */
2645    int numberSolutions = model->getSolutionCount();
2646    double cutoff = model->getCutoff();
2647    int method=0;
2648    int i;
2649    if (numberSolutions) {
2650      int numberHeuristic = model->getNumberHeuristicSolutions();
2651      if (numberHeuristic<numberSolutions) {
2652        method = 1;
2653      } else {
2654        method = 2;
2655        // look further
2656        for ( i = 0 ; i < numberObjects ; i++) {
2657          int numberNext = numberInfeasibilitiesUp[i];
2658         
2659          if (numberNext<numberUnsatisfied) {
2660            int numberUp = numberUnsatisfied - numberInfeasibilitiesUp[i];
2661            double perUnsatisfied = changeUp[i]/(double) numberUp;
2662            double estimatedObjective = objectiveValue + numberUnsatisfied * perUnsatisfied;
2663            if (estimatedObjective<cutoff) 
2664              method=3;
2665          }
2666          numberNext = numberInfeasibilitiesDown[i];
2667          if (numberNext<numberUnsatisfied) {
2668            int numberDown = numberUnsatisfied - numberInfeasibilitiesDown[i];
2669            double perUnsatisfied = changeDown[i]/(double) numberDown;
2670            double estimatedObjective = objectiveValue + numberUnsatisfied * perUnsatisfied;
2671            if (estimatedObjective<cutoff) 
2672              method=3;
2673          }
2674        }
2675      }
2676      method=2;
2677    } else {
2678      method = 0;
2679    }
2680    // Uncomment next to force method 4
2681    //method=4;
2682    /* Methods :
2683       0 - fewest infeasibilities
2684       1 - largest min change in objective
2685       2 - as 1 but use sum of changes if min close
2686       3 - predicted best solution
2687       4 - take cheapest up branch if infeasibilities same
2688    */
2689    int bestNumber=COIN_INT_MAX;
2690    double bestCriterion=-1.0e50;
2691    double alternativeCriterion = -1.0;
2692    double bestEstimate = 1.0e100;
2693    switch (method) {
2694    case 0:
2695      // could add in depth as well
2696      for ( i = 0 ; i < numberObjects ; i++) {
2697        int thisNumber = min(numberInfeasibilitiesUp[i],numberInfeasibilitiesDown[i]);
2698        if (thisNumber<=bestNumber) {
2699          int betterWay=0;
2700          if (numberInfeasibilitiesUp[i]<numberInfeasibilitiesDown[i]) {
2701            if (numberInfeasibilitiesUp[i]<bestNumber) {
2702              betterWay = 1;
2703            } else {
2704              if (changeUp[i]<bestCriterion)
2705                betterWay=1;
2706            }
2707          } else if (numberInfeasibilitiesUp[i]>numberInfeasibilitiesDown[i]) {
2708            if (numberInfeasibilitiesDown[i]<bestNumber) {
2709              betterWay = -1;
2710            } else {
2711              if (changeDown[i]<bestCriterion)
2712                betterWay=-1;
2713            }
2714          } else {
2715            // up and down have same number
2716            bool better=false;
2717            if (numberInfeasibilitiesUp[i]<bestNumber) {
2718              better=true;
2719            } else if (numberInfeasibilitiesUp[i]==bestNumber) {
2720              if (min(changeUp[i],changeDown[i])<bestCriterion)
2721                better=true;;
2722            }
2723            if (better) {
2724              // see which way
2725              if (changeUp[i]<=changeDown[i])
2726                betterWay=1;
2727              else
2728                betterWay=-1;
2729            }
2730          }
2731          if (betterWay) {
2732            bestCriterion = min(changeUp[i],changeDown[i]);
2733            bestNumber = thisNumber;
2734            whichObject = i;
2735            bestWay = betterWay;
2736          }
2737        }
2738      }
2739      break;
2740    case 1:
2741      for ( i = 0 ; i < numberObjects ; i++) {
2742        int betterWay=0;
2743        if (changeUp[i]<=changeDown[i]) {
2744          if (changeUp[i]>bestCriterion)
2745            betterWay=1;
2746        } else {
2747          if (changeDown[i]>bestCriterion)
2748            betterWay=-1;
2749        }
2750        if (betterWay) {
2751          bestCriterion = min(changeUp[i],changeDown[i]);
2752          whichObject = i;
2753          bestWay = betterWay;
2754        }
2755      }
2756      break;
2757    case 2:
2758      for ( i = 0 ; i < numberObjects ; i++) {
2759        double change = min(changeUp[i],changeDown[i]);
2760        double sum = changeUp[i] + changeDown[i];
2761        bool take=false;
2762        if (change>1.1*bestCriterion) 
2763          take=true;
2764        else if (change>0.9*bestCriterion&&sum+change>bestCriterion+alternativeCriterion) 
2765          take=true;
2766        if (take) {
2767          if (changeUp[i]<=changeDown[i]) {
2768            if (changeUp[i]>bestCriterion)
2769              bestWay=1;
2770          } else {
2771            if (changeDown[i]>bestCriterion)
2772              bestWay=-1;
2773          }
2774          bestCriterion = change;
2775          alternativeCriterion = sum;
2776          whichObject = i;
2777        }
2778      }
2779      break;
2780    case 3:
2781      for ( i = 0 ; i < numberObjects ; i++) {
2782        int numberNext = numberInfeasibilitiesUp[i];
2783       
2784        if (numberNext<numberUnsatisfied) {
2785          int numberUp = numberUnsatisfied - numberInfeasibilitiesUp[i];
2786          double perUnsatisfied = changeUp[i]/(double) numberUp;
2787          double estimatedObjective = objectiveValue + numberUnsatisfied * perUnsatisfied;
2788          if (estimatedObjective<bestEstimate) {
2789            bestEstimate = estimatedObjective;
2790            bestWay=1;
2791            whichObject=i;
2792          }
2793        }
2794        numberNext = numberInfeasibilitiesDown[i];
2795        if (numberNext<numberUnsatisfied) {
2796          int numberDown = numberUnsatisfied - numberInfeasibilitiesDown[i];
2797          double perUnsatisfied = changeDown[i]/(double) numberDown;
2798          double estimatedObjective = objectiveValue + numberUnsatisfied * perUnsatisfied;
2799          if (estimatedObjective<bestEstimate) {
2800            bestEstimate = estimatedObjective;
2801            bestWay=-1;
2802            whichObject=i;
2803          }
2804        }
2805      }
2806      break;
2807    case 4:
2808      // if number infeas same then cheapest up
2809      // first get best number or when going down
2810      // now choose smallest change up amongst equal number infeas
2811      for ( i = 0 ; i < numberObjects ; i++) {
2812        int thisNumber = min(numberInfeasibilitiesUp[i],numberInfeasibilitiesDown[i]);
2813        if (thisNumber<=bestNumber) {
2814          int betterWay=0;
2815          if (numberInfeasibilitiesUp[i]<numberInfeasibilitiesDown[i]) {
2816            if (numberInfeasibilitiesUp[i]<bestNumber) {
2817              betterWay = 1;
2818            } else {
2819              if (changeUp[i]<bestCriterion)
2820                betterWay=1;
2821            }
2822          } else if (numberInfeasibilitiesUp[i]>numberInfeasibilitiesDown[i]) {
2823            if (numberInfeasibilitiesDown[i]<bestNumber) {
2824              betterWay = -1;
2825            } else {
2826              if (changeDown[i]<bestCriterion)
2827                betterWay=-1;
2828            }
2829          } else {
2830            // up and down have same number
2831            bool better=false;
2832            if (numberInfeasibilitiesUp[i]<bestNumber) {
2833              better=true;
2834            } else if (numberInfeasibilitiesUp[i]==bestNumber) {
2835              if (min(changeUp[i],changeDown[i])<bestCriterion)
2836                better=true;;
2837            }
2838            if (better) {
2839              // see which way
2840              if (changeUp[i]<=changeDown[i])
2841                betterWay=1;
2842              else
2843                betterWay=-1;
2844            }
2845          }
2846          if (betterWay) {
2847            bestCriterion = min(changeUp[i],changeDown[i]);
2848            bestNumber = thisNumber;
2849            whichObject = i;
2850            bestWay = betterWay;
2851          }
2852        }
2853      }
2854      bestCriterion=1.0e50;
2855      for ( i = 0 ; i < numberObjects ; i++) {
2856        int thisNumber = numberInfeasibilitiesUp[i];
2857        if (thisNumber==bestNumber&&changeUp) {
2858          if (changeUp[i]<bestCriterion) {
2859            bestCriterion = changeUp[i];
2860            whichObject = i;
2861            bestWay = 1;
2862          }
2863        }
2864      }
2865      break;
2866    }
2867    // set way in best
2868    if (whichObject>=0) {
2869      CbcBranchingObject * bestObject = objects[whichObject];
2870      if (bestObject->object()&&bestObject->object()->preferredWay()) 
2871        bestWay = bestObject->object()->preferredWay();
2872      bestObject->way(bestWay);
2873    } else {
2874      printf("debug\n");
2875    }
2876  }
2877  return whichObject;
2878}
2879
2880//##############################################################################
2881
2882// Default Constructor
2883CbcFollowOn::CbcFollowOn ()
2884  : CbcObject(),
2885    rhs_(NULL)
2886{
2887}
2888
2889// Useful constructor
2890CbcFollowOn::CbcFollowOn (CbcModel * model)
2891  : CbcObject(model)
2892{
2893  assert (model);
2894  OsiSolverInterface * solver = model_->solver();
2895  matrix_ = *solver->getMatrixByCol();
2896  matrix_.removeGaps();
2897  matrix_.setExtraGap(0.0);
2898  matrixByRow_ = *solver->getMatrixByRow();
2899  int numberRows = matrix_.getNumRows();
2900 
2901  rhs_ = new int[numberRows];
2902  int i;
2903  const double * rowLower = solver->getRowLower();
2904  const double * rowUpper = solver->getRowUpper();
2905  // Row copy
2906  const double * elementByRow = matrixByRow_.getElements();
2907  const int * column = matrixByRow_.getIndices();
2908  const CoinBigIndex * rowStart = matrixByRow_.getVectorStarts();
2909  const int * rowLength = matrixByRow_.getVectorLengths();
2910  for (i=0;i<numberRows;i++) {
2911    rhs_[i]=0;
2912    double value = rowLower[i];
2913    if (value==rowUpper[i]) {
2914      if (floor(value)==value&&value>=1.0&&value<10.0) {
2915        // check elements
2916        bool good=true;
2917        for (int j=rowStart[i];j<rowStart[i]+rowLength[i];j++) {
2918          int iColumn = column[j];
2919          if (!solver->isBinary(iColumn))
2920            good=false;
2921          double elValue = elementByRow[j];
2922          if (floor(elValue)!=elValue||value<1.0)
2923            good=false;
2924        }
2925        if (good)
2926          rhs_[i]=(int) value;
2927      }
2928    }
2929  }
2930}
2931
2932// Copy constructor
2933CbcFollowOn::CbcFollowOn ( const CbcFollowOn & rhs)
2934  :CbcObject(rhs),
2935   matrix_(rhs.matrix_),
2936   matrixByRow_(rhs.matrixByRow_)
2937{
2938  int numberRows = matrix_.getNumRows();
2939  rhs_= CoinCopyOfArray(rhs.rhs_,numberRows);
2940}
2941
2942// Clone
2943CbcObject *
2944CbcFollowOn::clone() const
2945{
2946  return new CbcFollowOn(*this);
2947}
2948
2949// Assignment operator
2950CbcFollowOn & 
2951CbcFollowOn::operator=( const CbcFollowOn& rhs)
2952{
2953  if (this!=&rhs) {
2954    CbcObject::operator=(rhs);
2955    delete [] rhs_;
2956    matrix_ = rhs.matrix_;
2957    matrixByRow_ = rhs.matrixByRow_;
2958    int numberRows = matrix_.getNumRows();
2959    rhs_= CoinCopyOfArray(rhs.rhs_,numberRows);
2960  }
2961  return *this;
2962}
2963
2964// Destructor
2965CbcFollowOn::~CbcFollowOn ()
2966{
2967  delete [] rhs_;
2968}
2969// As some computation is needed in more than one place - returns row
2970int 
2971CbcFollowOn::gutsOfFollowOn(int & otherRow, int & preferredWay) const
2972{
2973  int whichRow=-1;
2974  otherRow=-1;
2975  int numberRows = matrix_.getNumRows();
2976 
2977  int i;
2978  // For sorting
2979  int * sort = new int [numberRows];
2980  int * isort = new int [numberRows];
2981  // Column copy
2982  //const double * element = matrix_.getElements();
2983  const int * row = matrix_.getIndices();
2984  const CoinBigIndex * columnStart = matrix_.getVectorStarts();
2985  const int * columnLength = matrix_.getVectorLengths();
2986  // Row copy
2987  const double * elementByRow = matrixByRow_.getElements();
2988  const int * column = matrixByRow_.getIndices();
2989  const CoinBigIndex * rowStart = matrixByRow_.getVectorStarts();
2990  const int * rowLength = matrixByRow_.getVectorLengths();
2991  OsiSolverInterface * solver = model_->solver();
2992  const double * columnLower = solver->getColLower();
2993  const double * columnUpper = solver->getColUpper();
2994  const double * solution = solver->getColSolution();
2995  double integerTolerance = model_->getDblParam(CbcModel::CbcIntegerTolerance);
2996  int nSort=0;
2997  for (i=0;i<numberRows;i++) {
2998    if (rhs_[i]) {
2999      // check elements
3000      double smallest=1.0e10;
3001      double largest=0.0;
3002      int rhsValue=rhs_[i];
3003      int number1=0;
3004      int numberUnsatisfied=0;
3005      for (int j=rowStart[i];j<rowStart[i]+rowLength[i];j++) {
3006        int iColumn = column[j];
3007        double value = elementByRow[j];
3008        double solValue = solution[iColumn];
3009        if (columnLower[iColumn]!=columnUpper[iColumn]) {
3010          smallest = CoinMin(smallest,value);
3011          largest = CoinMax(largest,value);
3012          if (value==1.0)
3013            number1++;
3014          if (solValue<1.0-integerTolerance&&solValue>integerTolerance)
3015            numberUnsatisfied++;
3016        } else {
3017          rhsValue -= (int)(value*floor(solValue+0.5));
3018        }
3019      }
3020      if (numberUnsatisfied>1) {
3021        if (smallest<largest) {
3022          // probably no good but check a few things
3023          assert (largest<=rhsValue);
3024          if (number1==1&&largest==rhsValue)
3025            printf("could fix\n");
3026        } else if (largest==rhsValue) {
3027          sort[nSort]=i;
3028          isort[nSort++]=-numberUnsatisfied;
3029        }
3030      }
3031    }
3032  }
3033  if (nSort>1) {
3034    CoinSort_2(isort,isort+nSort,sort);
3035    CoinZeroN(isort,numberRows);
3036    double * other = new double[numberRows];
3037    CoinZeroN(other,numberRows);
3038    int * which = new int[numberRows];
3039    //#define COUNT
3040#ifndef COUNT
3041    bool beforeSolution = model_->getSolutionCount()==0;
3042#endif
3043    for (int k=0;k<nSort-1;k++) {
3044      i=sort[k];
3045      int numberUnsatisfied = 0;
3046      int n=0;
3047      int j;
3048      for (j=rowStart[i];j<rowStart[i]+rowLength[i];j++) {
3049        int iColumn = column[j];
3050        if (columnLower[iColumn]!=columnUpper[iColumn]) {
3051          double solValue = solution[iColumn]-columnLower[iColumn];
3052          if (solValue<1.0-integerTolerance&&solValue>integerTolerance) {
3053            numberUnsatisfied++;
3054            for (int jj=columnStart[iColumn];jj<columnStart[iColumn]+columnLength[iColumn];jj++) {
3055              int iRow = row[jj];
3056              if (rhs_[iRow]) {
3057                other[iRow]+=solValue;
3058                if (isort[iRow]) {
3059                  isort[iRow]++;
3060                } else {
3061                  isort[iRow]=1;
3062                  which[n++]=iRow;
3063                }
3064              }
3065            }
3066          }
3067        }
3068      }
3069      double total=0.0;
3070      // Take out row
3071      double sumThis=other[i];
3072      other[i]=0.0;
3073      assert (numberUnsatisfied==isort[i]);
3074      // find one nearest half if solution, one if before solution
3075      int iBest=-1;
3076      double dtarget=0.5*total;
3077#ifdef COUNT
3078      int target = (numberUnsatisfied+1)>>1;
3079      int best=numberUnsatisfied;
3080#else
3081      double best;
3082      if (beforeSolution)
3083        best=dtarget;
3084      else
3085        best=1.0e30;
3086#endif
3087      for (j=0;j<n;j++) {
3088        int iRow = which[j];
3089        double dvalue=other[iRow];
3090        other[iRow]=0.0;
3091#ifdef COUNT
3092        int value = isort[iRow];
3093#endif
3094        isort[iRow]=0;
3095        if (fabs(dvalue)<1.0e-8||fabs(sumThis-dvalue)<1.0e-8)
3096          continue;
3097        if (dvalue<integerTolerance||dvalue>1.0-integerTolerance)
3098          continue;
3099#ifdef COUNT
3100        if (abs(value-target)<best&&value!=numberUnsatisfied) {
3101          best=abs(value-target);
3102          iBest=iRow;
3103          if (dvalue<dtarget)
3104            preferredWay=1;
3105          else
3106            preferredWay=-1;
3107        }
3108#else
3109        if (beforeSolution) {
3110          if (fabs(dvalue-dtarget)>best) {
3111            best = fabs(dvalue-dtarget);
3112            iBest=iRow;
3113            if (dvalue<dtarget)
3114              preferredWay=1;
3115            else
3116              preferredWay=-1;
3117          }
3118        } else {
3119          if (fabs(dvalue-dtarget)<best) {
3120            best = fabs(dvalue-dtarget);
3121            iBest=iRow;
3122            if (dvalue<dtarget)
3123              preferredWay=1;
3124            else
3125              preferredWay=-1;
3126          }
3127        }
3128#endif
3129      }
3130      if (iBest>=0) {
3131        whichRow=i;
3132        otherRow=iBest;
3133        break;
3134      }
3135    }
3136    delete [] which;
3137    delete [] other;
3138  }
3139  delete [] sort;
3140  delete [] isort;
3141  return whichRow;
3142}
3143
3144// Infeasibility - large is 0.5
3145double 
3146CbcFollowOn::infeasibility(int & preferredWay) const
3147{
3148  int otherRow=0;
3149  int whichRow = gutsOfFollowOn(otherRow,preferredWay);
3150  if (whichRow<0)
3151    return 0.0;
3152  else
3153  return 2.0* model_->getDblParam(CbcModel::CbcIntegerTolerance);
3154}
3155
3156// This looks at solution and sets bounds to contain solution
3157void 
3158CbcFollowOn::feasibleRegion()
3159{
3160}
3161
3162
3163// Creates a branching object
3164CbcBranchingObject * 
3165CbcFollowOn::createBranch(int way) 
3166{
3167  int otherRow=0;
3168  int preferredWay;
3169  int whichRow = gutsOfFollowOn(otherRow,preferredWay);
3170  assert(way==preferredWay);
3171  assert (whichRow>=0);
3172  int numberColumns = matrix_.getNumCols();
3173 
3174  // Column copy
3175  //const double * element = matrix_.getElements();
3176  const int * row = matrix_.getIndices();
3177  const CoinBigIndex * columnStart = matrix_.getVectorStarts();
3178  const int * columnLength = matrix_.getVectorLengths();
3179  // Row copy
3180  //const double * elementByRow = matrixByRow_.getElements();
3181  const int * column = matrixByRow_.getIndices();
3182  const CoinBigIndex * rowStart = matrixByRow_.getVectorStarts();
3183  const int * rowLength = matrixByRow_.getVectorLengths();
3184  OsiSolverInterface * solver = model_->solver();
3185  const double * columnLower = solver->getColLower();
3186  const double * columnUpper = solver->getColUpper();
3187  //const double * solution = solver->getColSolution();
3188  int nUp=0;
3189  int nDown=0;
3190  int * upList = new int[numberColumns];
3191  int * downList = new int[numberColumns];
3192  int j;
3193  for (j=rowStart[whichRow];j<rowStart[whichRow]+rowLength[whichRow];j++) {
3194    int iColumn = column[j];
3195    if (columnLower[iColumn]!=columnUpper[iColumn]) {
3196      bool up=true;
3197      for (int jj=columnStart[iColumn];jj<columnStart[iColumn]+columnLength[iColumn];jj++) {
3198        int iRow = row[jj];
3199        if (iRow==otherRow) {
3200          up=false;
3201          break;
3202        }
3203      }
3204      if (up)
3205        upList[nUp++]=iColumn;
3206      else
3207        downList[nDown++]=iColumn;
3208    }
3209  }
3210  //printf("way %d\n",way);
3211  // create object
3212  //printf("would fix %d down and %d up\n",nDown,nUp);
3213  CbcBranchingObject * branch
3214     = new CbcFixingBranchingObject(model_,way,
3215                                         nDown,downList,nUp,upList);
3216  delete [] upList;
3217  delete [] downList;
3218  return branch;
3219}
3220
3221//##############################################################################
3222
3223// Default Constructor
3224CbcFixingBranchingObject::CbcFixingBranchingObject()
3225  :CbcBranchingObject()
3226{
3227  numberDown_=0;
3228  numberUp_=0;
3229  downList_=NULL;
3230  upList_=NULL;
3231}
3232
3233// Useful constructor
3234CbcFixingBranchingObject::CbcFixingBranchingObject (CbcModel * model,
3235                                                    int way ,
3236                                                    int numberOnDownSide, const int * down,
3237                                                    int numberOnUpSide, const int * up)
3238  :CbcBranchingObject(model,0,way,0.5)
3239{
3240  numberDown_=numberOnDownSide;
3241  numberUp_=numberOnUpSide;
3242  downList_ = CoinCopyOfArray(down,numberDown_);
3243  upList_ = CoinCopyOfArray(up,numberUp_);
3244}
3245
3246// Copy constructor
3247CbcFixingBranchingObject::CbcFixingBranchingObject ( const CbcFixingBranchingObject & rhs) :CbcBranchingObject(rhs)
3248{
3249  numberDown_=rhs.numberDown_;
3250  numberUp_=rhs.numberUp_;
3251  downList_ = CoinCopyOfArray(rhs.downList_,numberDown_);
3252  upList_ = CoinCopyOfArray(rhs.upList_,numberUp_);
3253}
3254
3255// Assignment operator
3256CbcFixingBranchingObject & 
3257CbcFixingBranchingObject::operator=( const CbcFixingBranchingObject& rhs)
3258{
3259  if (this != &rhs) {
3260    CbcBranchingObject::operator=(rhs);
3261    delete [] downList_;
3262    delete [] upList_;
3263    numberDown_=rhs.numberDown_;
3264    numberUp_=rhs.numberUp_;
3265    downList_ = CoinCopyOfArray(rhs.downList_,numberDown_);
3266    upList_ = CoinCopyOfArray(rhs.upList_,numberUp_);
3267  }
3268  return *this;
3269}
3270CbcBranchingObject * 
3271CbcFixingBranchingObject::clone() const
3272{ 
3273  return (new CbcFixingBranchingObject(*this));
3274}
3275
3276
3277// Destructor
3278CbcFixingBranchingObject::~CbcFixingBranchingObject ()
3279{
3280  delete [] downList_;
3281  delete [] upList_;
3282}
3283double
3284CbcFixingBranchingObject::branch()
3285{
3286  decrementNumberBranchesLeft();
3287  OsiSolverInterface * solver = model_->solver();
3288  const double * columnLower = solver->getColLower();
3289  int i;
3290  // *** for way - up means fix all those in up section
3291  if (way_<0) {
3292#ifdef FULL_PRINT
3293    printf("Down Fix ");
3294#endif
3295    //printf("Down Fix %d\n",numberDown_);
3296    for (i=0;i<numberDown_;i++) {
3297      int iColumn = downList_[i];
3298      model_->solver()->setColUpper(iColumn,columnLower[iColumn]);
3299#ifdef FULL_PRINT
3300      printf("Setting bound on %d to lower bound\n",iColumn);
3301#endif
3302    }
3303    way_=1;       // Swap direction
3304  } else {
3305#ifdef FULL_PRINT
3306    printf("Up Fix ");
3307#endif
3308    //printf("Up Fix %d\n",numberUp_);
3309    for (i=0;i<numberUp_;i++) {
3310      int iColumn = upList_[i];
3311      model_->solver()->setColUpper(iColumn,columnLower[iColumn]);
3312#ifdef FULL_PRINT
3313      printf("Setting bound on %d to lower bound\n",iColumn);
3314#endif
3315    }
3316    way_=-1;      // Swap direction
3317  }
3318#ifdef FULL_PRINT
3319  printf("\n");
3320#endif
3321  return 0.0;
3322}
3323void
3324CbcFixingBranchingObject::print()
3325{
3326  int i;
3327  // *** for way - up means fix all those in up section
3328  if (way_<0) {
3329    printf("Down Fix ");
3330    for (i=0;i<numberDown_;i++) {
3331      int iColumn = downList_[i];
3332      printf("%d ",iColumn);
3333    }
3334  } else {
3335    printf("Up Fix ");
3336    for (i=0;i<numberUp_;i++) {
3337      int iColumn = upList_[i];
3338      printf("%d ",iColumn);
3339    }
3340  }
3341  printf("\n");
3342}
3343
3344/** Compare the original object of \c this with the original object of \c
3345    brObj. Assumes that there is an ordering of the original objects.
3346    This method should be invoked only if \c this and brObj are of the same
3347    type.
3348    Return negative/0/positive depending on whether \c this is
3349    smaller/same/larger than the argument.
3350*/
3351int
3352CbcFixingBranchingObject::compareOriginalObject
3353(const CbcBranchingObject* brObj) const
3354{
3355  throw("must implement");
3356}
3357
3358/** Compare the \c this with \c brObj. \c this and \c brObj must be os the
3359    same type and must have the same original object, but they may have
3360    different feasible regions.
3361    Return the appropriate CbcRangeCompare value (first argument being the
3362    sub/superset if that's the case). In case of overlap (and if \c
3363    replaceIfOverlap is true) replace the current branching object with one
3364    whose feasible region is the overlap.
3365   */
3366CbcRangeCompare
3367CbcFixingBranchingObject::compareBranchingObject
3368(const CbcBranchingObject* brObj, const bool replaceIfOverlap)
3369{
3370  const CbcFixingBranchingObject* br =
3371    dynamic_cast<const CbcFixingBranchingObject*>(brObj);
3372  assert(br);
3373  // If two FixingBranchingObject's have the same base object then it's pretty
3374  // much guaranteed
3375  throw("must implement");
3376}
3377
3378//##############################################################################
3379
3380// Default Constructor
3381CbcNWay::CbcNWay ()
3382  : CbcObject(),
3383    numberMembers_(0),
3384    members_(NULL),
3385    consequence_(NULL)
3386{
3387}
3388
3389// Useful constructor (which are integer indices)
3390CbcNWay::CbcNWay (CbcModel * model, int numberMembers,
3391                  const int * which, int identifier)
3392  : CbcObject(model)
3393{
3394  id_=identifier;
3395  numberMembers_=numberMembers;
3396  if (numberMembers_) {
3397    members_ = new int[numberMembers_];
3398    memcpy(members_,which,numberMembers_*sizeof(int));
3399  } else {
3400    members_ = NULL;
3401  }
3402  consequence_ = NULL;
3403}
3404
3405// Copy constructor
3406CbcNWay::CbcNWay ( const CbcNWay & rhs)
3407  :CbcObject(rhs)
3408{
3409  numberMembers_ = rhs.numberMembers_;
3410  consequence_ = NULL;
3411  if (numberMembers_) {
3412    members_ = new int[numberMembers_];
3413    memcpy(members_,rhs.members_,numberMembers_*sizeof(int));
3414    if (rhs.consequence_) {
3415      consequence_ = new CbcConsequence * [numberMembers_];
3416      for (int i=0;i<numberMembers_;i++) {
3417        if (rhs.consequence_[i])
3418          consequence_[i]= rhs.consequence_[i]->clone();
3419        else
3420          consequence_[i]=NULL;
3421      }
3422    }
3423  } else {
3424    members_ = NULL;
3425  }
3426}
3427
3428// Clone
3429CbcObject *
3430CbcNWay::clone() const
3431{
3432  return new CbcNWay(*this);
3433}
3434
3435// Assignment operator
3436CbcNWay & 
3437CbcNWay::operator=( const CbcNWay& rhs)
3438{
3439  if (this!=&rhs) {
3440    CbcObject::operator=(rhs);
3441    delete [] members_;
3442    numberMembers_ = rhs.numberMembers_;
3443    if (consequence_) {
3444      for (int i=0;i<numberMembers_;i++) 
3445        delete consequence_[i];
3446      delete [] consequence_;
3447      consequence_=NULL;
3448    }
3449    if (numberMembers_) {
3450      members_ = new int[numberMembers_];
3451      memcpy(members_,rhs.members_,numberMembers_*sizeof(int));
3452    } else {
3453      members_ = NULL;
3454    }
3455    if (rhs.consequence_) {
3456      consequence_ = new CbcConsequence * [numberMembers_];
3457      for (int i=0;i<numberMembers_;i++) {
3458        if (rhs.consequence_[i])
3459          consequence_[i]= rhs.consequence_[i]->clone();
3460        else
3461          consequence_[i]=NULL;
3462      }
3463    }
3464  }
3465  return *this;
3466}
3467
3468// Destructor
3469CbcNWay::~CbcNWay ()
3470{
3471  delete [] members_;
3472  if (consequence_) {
3473    for (int i=0;i<numberMembers_;i++) 
3474      delete consequence_[i];
3475    delete [] consequence_;
3476  }
3477}
3478// Set up a consequence for a single member
3479void 
3480CbcNWay::setConsequence(int iColumn, const CbcConsequence & consequence)
3481{
3482  if (!consequence_) {
3483    consequence_ = new CbcConsequence * [numberMembers_];
3484    for (int i=0;i<numberMembers_;i++) 
3485      consequence_[i]=NULL;
3486  }
3487  for (int i=0;i<numberMembers_;i++) {
3488    if (members_[i]==iColumn) {
3489      consequence_[i]=consequence.clone();
3490      break;
3491    }
3492  }
3493}
3494
3495// Applies a consequence for a single member
3496void 
3497CbcNWay::applyConsequence(int iSequence, int state) const
3498{
3499  assert (state==-9999||state==9999);
3500  if (consequence_) {
3501    CbcConsequence * consequence = consequence_[iSequence];
3502    if (consequence) 
3503      consequence->applyToSolver(model_->solver(),state);
3504  }
3505}
3506 
3507// Infeasibility - large is 0.5
3508double 
3509CbcNWay::infeasibility(int & preferredWay) const
3510{
3511  int numberUnsatis=0;
3512  int j;
3513  OsiSolverInterface * solver = model_->solver();
3514  const double * solution = model_->testSolution();
3515  const double * lower = solver->getColLower();
3516  const double * upper = solver->getColUpper();
3517  double largestValue=0.0;
3518 
3519  double integerTolerance = 
3520    model_->getDblParam(CbcModel::CbcIntegerTolerance);
3521
3522  for (j=0;j<numberMembers_;j++) {
3523    int iColumn = members_[j];
3524    double value = solution[iColumn];
3525    value = CoinMax(value, lower[iColumn]);
3526    value = CoinMin(value, upper[iColumn]);
3527    double distance = CoinMin(value-lower[iColumn],upper[iColumn]-value);
3528    if (distance>integerTolerance) {
3529      numberUnsatis++;
3530      largestValue = CoinMax(distance,largestValue);
3531    }
3532  }
3533  preferredWay=1;
3534  if (numberUnsatis) {
3535    return largestValue;
3536  } else {
3537    return 0.0; // satisfied
3538  }
3539}
3540
3541// This looks at solution and sets bounds to contain solution
3542void 
3543CbcNWay::feasibleRegion()
3544{
3545  int j;
3546  OsiSolverInterface * solver = model_->solver();
3547  const double * solution = model_->testSolution();
3548  const double * lower = solver->getColLower();
3549  const double * upper = solver->getColUpper();
3550  double integerTolerance = 
3551    model_->getDblParam(CbcModel::CbcIntegerTolerance);
3552  for (j=0;j<numberMembers_;j++) {
3553    int iColumn = members_[j];
3554    double value = solution[iColumn];
3555    value = CoinMax(value, lower[iColumn]);
3556    value = CoinMin(value, upper[iColumn]);
3557    if (value>=upper[iColumn]-integerTolerance) {
3558      solver->setColLower(iColumn,upper[iColumn]);
3559    } else {
3560      assert (value<=lower[iColumn]+integerTolerance);
3561      solver->setColUpper(iColumn,lower[iColumn]);
3562    }
3563  }
3564}
3565// Redoes data when sequence numbers change
3566void 
3567CbcNWay::redoSequenceEtc(CbcModel * model, int numberColumns, const int * originalColumns)
3568{
3569  model_=model;
3570  int n2=0;
3571  for (int j=0;j<numberMembers_;j++) {
3572    int iColumn = members_[j];
3573    int i;
3574    for (i=0;i<numberColumns;i++) {
3575      if (originalColumns[i]==iColumn)
3576        break;
3577    }
3578    if (i<numberColumns) {
3579      members_[n2]=i;
3580      consequence_[n2++]=consequence_[j];
3581    } else {
3582      delete consequence_[j];
3583    }
3584  }
3585  if (n2<numberMembers_) {
3586    printf("** NWay number of members reduced from %d to %d!\n",numberMembers_,n2);
3587    numberMembers_=n2;
3588  }
3589}
3590
3591
3592// Creates a branching object
3593CbcBranchingObject * 
3594CbcNWay::createBranch(int way) 
3595{
3596  int numberFree=0;
3597  int j;
3598
3599  OsiSolverInterface * solver = model_->solver();
3600  const double * solution = model_->testSolution();
3601  const double * lower = solver->getColLower();
3602  const double * upper = solver->getColUpper();
3603  int * list = new int[numberMembers_];
3604  double * sort = new double[numberMembers_];
3605
3606  for (j=0;j<numberMembers_;j++) {
3607    int iColumn = members_[j];
3608    double value = solution[iColumn];
3609    value = CoinMax(value, lower[iColumn]);
3610    value = CoinMin(value, upper[iColumn]);
3611    if (upper[iColumn]>lower[iColumn]) {
3612      double distance = upper[iColumn]-value;
3613      list[numberFree]=j;
3614      sort[numberFree++]=distance;
3615    }
3616  }
3617  assert (numberFree);
3618  // sort
3619  CoinSort_2(sort,sort+numberFree,list);
3620  // create object
3621  CbcBranchingObject * branch;
3622  branch = new CbcNWayBranchingObject(model_,this,numberFree,list);
3623  branch->setOriginalObject(this);
3624  delete [] list;
3625  delete [] sort;
3626  return branch;
3627}
3628 
3629//##############################################################################
3630
3631// Default Constructor
3632CbcNWayBranchingObject::CbcNWayBranchingObject()
3633  :CbcBranchingObject()
3634{
3635  order_=NULL;
3636  object_=NULL;
3637  numberInSet_=0;
3638  way_=0;
3639}
3640
3641// Useful constructor
3642CbcNWayBranchingObject::CbcNWayBranchingObject (CbcModel * model,
3643                                                const CbcNWay * nway, 
3644                                                int number, const int * order)
3645  :CbcBranchingObject(model,nway->id(),-1,0.5)
3646{
3647  numberBranches_ = number;
3648  order_ = new int [number];
3649  object_=nway;
3650  numberInSet_=number;
3651  memcpy(order_,order,number*sizeof(int));
3652}
3653
3654// Copy constructor
3655CbcNWayBranchingObject::CbcNWayBranchingObject ( const CbcNWayBranchingObject & rhs) :CbcBranchingObject(rhs)
3656{
3657  numberInSet_=rhs.numberInSet_;
3658  object_=rhs.object_;
3659  if (numberInSet_) {
3660    order_ = new int [numberInSet_];
3661    memcpy(order_,rhs.order_,numberInSet_*sizeof(int));
3662  } else {
3663    order_=NULL;
3664  }   
3665}
3666
3667// Assignment operator
3668CbcNWayBranchingObject & 
3669CbcNWayBranchingObject::operator=( const CbcNWayBranchingObject& rhs)
3670{
3671  if (this != &rhs) {
3672    CbcBranchingObject::operator=(rhs);
3673    object_=rhs.object_;
3674    delete [] order_;
3675    numberInSet_=rhs.numberInSet_;
3676    if (numberInSet_) {
3677      order_ = new int [numberInSet_];
3678      memcpy(order_,rhs.order_,numberInSet_*sizeof(int));
3679    } else {
3680      order_=NULL;
3681    }   
3682  }
3683  return *this;
3684}
3685CbcBranchingObject * 
3686CbcNWayBranchingObject::clone() const
3687{ 
3688  return (new CbcNWayBranchingObject(*this));
3689}
3690
3691
3692// Destructor
3693CbcNWayBranchingObject::~CbcNWayBranchingObject ()
3694{
3695  delete [] order_;
3696}
3697double
3698CbcNWayBranchingObject::branch()
3699{
3700  int which = branchIndex_;
3701  branchIndex_++;
3702  assert (numberBranchesLeft()>=0);
3703  if (which==0) {
3704    // first branch so way_ may mean something
3705    assert (way_==-1||way_==1);
3706    if (way_==-1)
3707      which++;
3708  } else if (which==1) {
3709    // second branch so way_ may mean something
3710    assert (way_==-1||way_==1);
3711    if (way_==-1)
3712      which--;
3713    // switch way off
3714    way_=0;
3715  }
3716  const double * lower = model_->solver()->getColLower();
3717  const double * upper = model_->solver()->getColUpper();
3718  const int * members = object_->members();
3719  for (int j=0;j<numberInSet_;j++) {
3720    int iSequence = order_[j];
3721    int iColumn = members[iSequence];
3722    if (j!=which) {
3723      model_->solver()->setColUpper(iColumn,lower[iColumn]);
3724      //model_->solver()->setColLower(iColumn,lower[iColumn]);
3725      assert (lower[iColumn]>-1.0e20);
3726      // apply any consequences
3727      object_->applyConsequence(iSequence,-9999);
3728    } else {
3729      model_->solver()->setColLower(iColumn,upper[iColumn]);
3730      //model_->solver()->setColUpper(iColumn,upper[iColumn]);
3731#ifdef FULL_PRINT
3732      printf("Up Fix %d to %g\n",iColumn,upper[iColumn]);
3733#endif
3734      assert (upper[iColumn]<1.0e20);
3735      // apply any consequences
3736      object_->applyConsequence(iSequence,9999);
3737    }
3738  }
3739  return 0.0;
3740}
3741void
3742CbcNWayBranchingObject::print()
3743{
3744  printf("NWay - Up Fix ");
3745  const int * members = object_->members();
3746  for (int j=0;j<way_;j++) {
3747    int iColumn = members[order_[j]];
3748    printf("%d ",iColumn);
3749  }
3750  printf("\n");
3751}
3752
3753/** Compare the original object of \c this with the original object of \c
3754    brObj. Assumes that there is an ordering of the original objects.
3755    This method should be invoked only if \c this and brObj are of the same
3756    type.
3757    Return negative/0/positive depending on whether \c this is
3758    smaller/same/larger than the argument.
3759*/
3760int
3761CbcNWayBranchingObject::compareOriginalObject
3762(const CbcBranchingObject* brObj) const
3763{
3764  throw("must implement");
3765}
3766
3767/** Compare the \c this with \c brObj. \c this and \c brObj must be os the
3768    same type and must have the same original object, but they may have
3769    different feasible regions.
3770    Return the appropriate CbcRangeCompare value (first argument being the
3771    sub/superset if that's the case). In case of overlap (and if \c
3772    replaceIfOverlap is true) replace the current branching object with one
3773    whose feasible region is the overlap.
3774*/
3775CbcRangeCompare
3776CbcNWayBranchingObject::compareBranchingObject
3777(const CbcBranchingObject* brObj, const bool replaceIfOverlap)
3778{
3779  throw("must implement");
3780}
3781
3782//##############################################################################
3783
3784// Default Constructor
3785CbcFixVariable::CbcFixVariable ()
3786  : CbcConsequence(),
3787    numberStates_(0),
3788    states_(NULL),
3789    startLower_(NULL),
3790    startUpper_(NULL),
3791    newBound_(NULL),
3792    variable_(NULL)
3793{
3794}
3795
3796// One useful Constructor
3797CbcFixVariable::CbcFixVariable (int numberStates,const int * states, const int * numberNewLower, 
3798                                const int ** newLowerValue,
3799                                const int ** lowerColumn,
3800                                const int * numberNewUpper, const int ** newUpperValue,
3801                                const int ** upperColumn)
3802  : CbcConsequence(),
3803    states_(NULL),
3804    startLower_(NULL),
3805    startUpper_(NULL),
3806    newBound_(NULL),
3807    variable_(NULL)
3808{
3809  // How much space
3810  numberStates_ = numberStates;
3811  if (numberStates_) {
3812    states_ = new int[numberStates_];
3813    memcpy(states_,states,numberStates_*sizeof(int));
3814    int i;
3815    int n=0;
3816    startLower_ = new int[numberStates_+1];
3817    startUpper_ = new int[numberStates_+1];
3818    startLower_[0]=0;
3819    //count
3820    for (i=0;i<numberStates_;i++) {
3821      n += numberNewLower[i];
3822      startUpper_[i]=n;
3823      n += numberNewUpper[i];
3824      startLower_[i+1]=n;
3825    }
3826    newBound_ = new double [n];
3827    variable_ = new int [n];
3828    n=0;
3829    for (i=0;i<numberStates_;i++) {
3830      int j;
3831      int k;
3832      const int * bound;
3833      const int * variable;
3834      k=numberNewLower[i];
3835      bound = newLowerValue[i];
3836      variable = lowerColumn[i];
3837      for (j=0;j<k;j++) {
3838        newBound_[n]=bound[j];
3839        variable_[n++]=variable[j];
3840      }
3841      k=numberNewUpper[i];
3842      bound = newUpperValue[i];
3843      variable = upperColumn[i];
3844      for (j=0;j<k;j++) {
3845        newBound_[n]=bound[j];
3846        variable_[n++]=variable[j];
3847      }
3848    }
3849  }
3850}
3851
3852// Copy constructor
3853CbcFixVariable::CbcFixVariable ( const CbcFixVariable & rhs)
3854  :CbcConsequence(rhs)
3855{
3856  numberStates_ = rhs.numberStates_;
3857  states_ = NULL;
3858  startLower_ = NULL;
3859  startUpper_ = NULL;
3860  newBound_ = NULL;
3861  variable_ = NULL;
3862  if (numberStates_) {
3863    states_ = CoinCopyOfArray(rhs.states_,numberStates_);
3864    startLower_ = CoinCopyOfArray(rhs.startLower_,numberStates_+1);
3865    startUpper_ = CoinCopyOfArray(rhs.startUpper_,numberStates_+1);
3866    int n=startLower_[numberStates_];
3867    newBound_ = CoinCopyOfArray(rhs.newBound_,n);
3868    variable_ = CoinCopyOfArray(rhs.variable_,n);
3869  }
3870}
3871
3872// Clone
3873CbcConsequence *
3874CbcFixVariable::clone() const
3875{
3876  return new CbcFixVariable(*this);
3877}
3878
3879// Assignment operator
3880CbcFixVariable & 
3881CbcFixVariable::operator=( const CbcFixVariable& rhs)
3882{
3883  if (this!=&rhs) {
3884    CbcConsequence::operator=(rhs);
3885    delete [] states_;
3886    delete [] startLower_;
3887    delete [] startUpper_;
3888    delete [] newBound_;
3889    delete [] variable_;
3890    states_ = NULL;
3891    startLower_ = NULL;
3892    startUpper_ = NULL;
3893    newBound_ = NULL;
3894    variable_ = NULL;
3895    numberStates_ = rhs.numberStates_;
3896    if (numberStates_) {
3897      states_ = CoinCopyOfArray(rhs.states_,numberStates_);
3898      startLower_ = CoinCopyOfArray(rhs.startLower_,numberStates_+1);
3899      startUpper_ = CoinCopyOfArray(rhs.startUpper_,numberStates_+1);
3900      int n=startLower_[numberStates_];
3901      newBound_ = CoinCopyOfArray(rhs.newBound_,n);
3902      variable_ = CoinCopyOfArray(rhs.variable_,n);
3903    }
3904  }
3905  return *this;
3906}
3907
3908// Destructor
3909CbcFixVariable::~CbcFixVariable ()
3910{
3911  delete [] states_;
3912  delete [] startLower_;
3913  delete [] startUpper_;
3914  delete [] newBound_;
3915  delete [] variable_;
3916}
3917// Set up a startLower for a single member
3918void 
3919CbcFixVariable::applyToSolver(OsiSolverInterface * solver, int state) const
3920{
3921  assert (state==-9999||state==9999);
3922  // Find state
3923  int find;
3924  for (find=0;find<numberStates_;find++) 
3925    if (states_[find]==state)
3926      break;
3927  if (find==numberStates_)
3928    return;
3929  int i;
3930  // Set new lower bounds
3931  for (i=startLower_[find];i<startUpper_[find];i++) {
3932    int iColumn = variable_[i];
3933    double value = newBound_[i];
3934    double oldValue = solver->getColLower()[iColumn];
3935    //printf("for %d old lower bound %g, new %g",iColumn,oldValue,value);
3936    solver->setColLower(iColumn,CoinMax(value,oldValue));
3937    //printf(" => %g\n",solver->getColLower()[iColumn]);
3938  }
3939  // Set new upper bounds
3940  for (i=startUpper_[find];i<startLower_[find+1];i++) {
3941    int iColumn = variable_[i];
3942    double value = newBound_[i];
3943    double oldValue = solver->getColUpper()[iColumn];
3944    //printf("for %d old upper bound %g, new %g",iColumn,oldValue,value);
3945    solver->setColUpper(iColumn,CoinMin(value,oldValue));
3946    //printf(" => %g\n",solver->getColUpper()[iColumn]);
3947  }
3948}
3949
3950//##############################################################################
3951
3952// Default Constructor
3953CbcDummyBranchingObject::CbcDummyBranchingObject(CbcModel * model)
3954  :CbcBranchingObject(model,0,0,0.5)
3955{
3956  setNumberBranchesLeft(1);
3957}
3958
3959
3960// Copy constructor
3961CbcDummyBranchingObject::CbcDummyBranchingObject ( const CbcDummyBranchingObject & rhs) :CbcBranchingObject(rhs)
3962{
3963}
3964
3965// Assignment operator
3966CbcDummyBranchingObject & 
3967CbcDummyBranchingObject::operator=( const CbcDummyBranchingObject& rhs)
3968{
3969  if (this != &rhs) {
3970    CbcBranchingObject::operator=(rhs);
3971  }
3972  return *this;
3973}
3974CbcBranchingObject * 
3975CbcDummyBranchingObject::clone() const
3976{ 
3977  return (new CbcDummyBranchingObject(*this));
3978}
3979
3980
3981// Destructor
3982CbcDummyBranchingObject::~CbcDummyBranchingObject ()
3983{
3984}
3985
3986/*
3987  Perform a dummy branch
3988*/
3989double
3990CbcDummyBranchingObject::branch()
3991{
3992  decrementNumberBranchesLeft();
3993  return 0.0;
3994}
3995// Print what would happen 
3996void
3997CbcDummyBranchingObject::print()
3998{
3999  printf("Dummy branch\n");
4000}
4001
4002/** Compare the original object of \c this with the original object of \c
4003    brObj. Assumes that there is an ordering of the original objects.
4004    This method should be invoked only if \c this and brObj are of the same
4005    type.
4006    Return negative/0/positive depending on whether \c this is
4007    smaller/same/larger than the argument.
4008*/
4009int
4010CbcDummyBranchingObject::compareOriginalObject
4011(const CbcBranchingObject* brObj) const
4012{
4013  throw("must implement");
4014}
4015
4016/** Compare the \c this with \c brObj. \c this and \c brObj must be os the
4017    same type and must have the same original object, but they may have
4018    different feasible regions.
4019    Return the appropriate CbcRangeCompare value (first argument being the
4020    sub/superset if that's the case). In case of overlap (and if \c
4021    replaceIfOverlap is true) replace the current branching object with one
4022    whose feasible region is the overlap.
4023*/
4024CbcRangeCompare
4025CbcDummyBranchingObject::compareBranchingObject
4026(const CbcBranchingObject* brObj, const bool replaceIfOverlap)
4027{
4028  throw("must implement");
4029}
4030
Note: See TracBrowser for help on using the repository browser.