source: branches/devel/Cbc/src/CbcBranchActual.cpp @ 417

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

for allCuts

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