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

Last change on this file since 720 was 720, checked in by forrest, 13 years ago

message handling and sos to CbcSolver?

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