source: trunk/CbcBranchActual.cpp @ 216

Last change on this file since 216 was 216, checked in by forrest, 15 years ago

add branching stuff

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