source: trunk/CbcBranchDynamic.cpp @ 307

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

for bonmin

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 27.7 KB
Line 
1// Copyright (C) 2002, International Business Machines
2// Corporation and others.  All Rights Reserved.
3#if defined(_MSC_VER)
4// Turn off compiler warning about long names
5#  pragma warning(disable:4786)
6#endif
7#include <cassert>
8#include <cmath>
9#include <cfloat>
10//#define CBC_DEBUG
11
12#include "OsiSolverInterface.hpp"
13#include "OsiSolverBranch.hpp"
14#include "CbcModel.hpp"
15#include "CbcMessage.hpp"
16#include "CbcBranchDynamic.hpp"
17#include "CoinSort.hpp"
18#include "CoinError.hpp"
19
20/** Default Constructor
21
22  Equivalent to an unspecified binary variable.
23*/
24CbcSimpleIntegerDynamicPseudoCost::CbcSimpleIntegerDynamicPseudoCost ()
25  : CbcSimpleInteger(),
26    downDynamicPseudoCost_(1.0e-5),
27    upDynamicPseudoCost_(1.0e-5),
28    upDownSeparator_(-1.0),
29    sumDownCost_(0.0),
30    sumUpCost_(0.0),
31    sumDownChange_(0.0),
32    sumUpChange_(0.0),
33    sumDownCostSquared_(0.0),
34    sumUpCostSquared_(0.0),
35    sumDownDecrease_(0.0),
36    sumUpDecrease_(0.0),
37    lastDownCost_(0.0),
38    lastUpCost_(0.0),
39    lastDownDecrease_(0),
40    lastUpDecrease_(0),
41    numberTimesDown_(0),
42    numberTimesUp_(0),
43    numberTimesDownInfeasible_(0),
44    numberTimesUpInfeasible_(0),
45    numberBeforeTrust_(0),
46    method_(0)
47{
48}
49
50/** Useful constructor
51
52  Loads dynamic upper & lower bounds for the specified variable.
53*/
54CbcSimpleIntegerDynamicPseudoCost::CbcSimpleIntegerDynamicPseudoCost (CbcModel * model, int sequence,
55                                    int iColumn, double breakEven)
56  : CbcSimpleInteger(model,sequence,iColumn,breakEven),
57    upDownSeparator_(-1.0),
58    sumDownCost_(0.0),
59    sumUpCost_(0.0),
60    sumDownChange_(0.0),
61    sumUpChange_(0.0),
62    sumDownCostSquared_(0.0),
63    sumUpCostSquared_(0.0),
64    sumDownDecrease_(0.0),
65    sumUpDecrease_(0.0),
66    lastDownCost_(0.0),
67    lastUpCost_(0.0),
68    lastDownDecrease_(0),
69    lastUpDecrease_(0),
70    numberTimesDown_(0),
71    numberTimesUp_(0),
72    numberTimesDownInfeasible_(0),
73    numberTimesUpInfeasible_(0),
74    numberBeforeTrust_(0),
75    method_(0)
76{
77  const double * cost = model->getObjCoefficients();
78  double costValue = CoinMax(1.0e-5,fabs(cost[iColumn]));
79  // treat as if will cost what it says up
80  upDynamicPseudoCost_=costValue;
81  // and balance at breakeven
82  downDynamicPseudoCost_=((1.0-breakEven_)*upDynamicPseudoCost_)/breakEven_;
83  // so initial will have some effect
84  sumUpCost_ = 2.0*upDynamicPseudoCost_;
85  sumUpChange_ = 2.0;
86  numberTimesUp_ = 2;
87  sumDownCost_ = 2.0*downDynamicPseudoCost_;
88  sumDownChange_ = 2.0;
89  numberTimesDown_ = 2;
90#if 0
91  // No
92  sumUpCost_ = 0.0;
93  sumUpChange_ = 0.0;
94  numberTimesUp_ = 0;
95  sumDownCost_ = 0.0;
96  sumDownChange_ = 0.0;
97  numberTimesDown_ = 0;
98#else
99  sumUpCost_ = 1.0*upDynamicPseudoCost_;
100  sumUpChange_ = 1.0;
101  numberTimesUp_ = 1;
102  sumDownCost_ = 1.0*downDynamicPseudoCost_;
103  sumDownChange_ = 1.0;
104  numberTimesDown_ = 1;
105#endif
106}
107
108/** Useful constructor
109
110  Loads dynamic upper & lower bounds for the specified variable.
111*/
112CbcSimpleIntegerDynamicPseudoCost::CbcSimpleIntegerDynamicPseudoCost (CbcModel * model, int sequence,
113                                    int iColumn, double downDynamicPseudoCost,
114                                                        double upDynamicPseudoCost)
115  : CbcSimpleInteger(model,sequence,iColumn),
116    upDownSeparator_(-1.0),
117    sumDownCost_(0.0),
118    sumUpCost_(0.0),
119    sumDownChange_(0.0),
120    sumUpChange_(0.0),
121    sumDownCostSquared_(0.0),
122    sumUpCostSquared_(0.0),
123    sumDownDecrease_(0.0),
124    sumUpDecrease_(0.0),
125    lastDownCost_(0.0),
126    lastUpCost_(0.0),
127    lastDownDecrease_(0),
128    lastUpDecrease_(0),
129    numberTimesDown_(0),
130    numberTimesUp_(0),
131    numberTimesDownInfeasible_(0),
132    numberTimesUpInfeasible_(0),
133    numberBeforeTrust_(0),
134    method_(0)
135{
136  downDynamicPseudoCost_ = downDynamicPseudoCost;
137  upDynamicPseudoCost_ = upDynamicPseudoCost;
138  breakEven_ = upDynamicPseudoCost_/(upDynamicPseudoCost_+downDynamicPseudoCost_);
139  // so initial will have some effect
140  sumUpCost_ = 2.0*upDynamicPseudoCost_;
141  sumUpChange_ = 2.0;
142  numberTimesUp_ = 2;
143  sumDownCost_ = 2.0*downDynamicPseudoCost_;
144  sumDownChange_ = 2.0;
145  numberTimesDown_ = 2;
146#if 0
147  // No
148  sumUpCost_ = 0.0;
149  sumUpChange_ = 0.0;
150  numberTimesUp_ = 0;
151  sumDownCost_ = 0.0;
152  sumDownChange_ = 0.0;
153  numberTimesDown_ = 0;
154#else
155  sumUpCost_ = 1.0*upDynamicPseudoCost_;
156  sumUpChange_ = 1.0;
157  numberTimesUp_ = 1;
158  sumDownCost_ = 1.0*downDynamicPseudoCost_;
159  sumDownChange_ = 1.0;
160  numberTimesDown_ = 1;
161#endif
162}
163
164// Copy constructor
165CbcSimpleIntegerDynamicPseudoCost::CbcSimpleIntegerDynamicPseudoCost ( const CbcSimpleIntegerDynamicPseudoCost & rhs)
166  :CbcSimpleInteger(rhs),
167   downDynamicPseudoCost_(rhs.downDynamicPseudoCost_),
168   upDynamicPseudoCost_(rhs.upDynamicPseudoCost_),
169   upDownSeparator_(rhs.upDownSeparator_),
170   sumDownCost_(rhs.sumDownCost_),
171   sumUpCost_(rhs.sumUpCost_),
172   sumDownChange_(rhs.sumDownChange_),
173   sumUpChange_(rhs.sumUpChange_),
174   sumDownCostSquared_(rhs.sumDownCostSquared_),
175   sumUpCostSquared_(rhs.sumUpCostSquared_),
176   sumDownDecrease_(rhs.sumDownDecrease_),
177   sumUpDecrease_(rhs.sumUpDecrease_),
178   lastDownCost_(rhs.lastDownCost_),
179   lastUpCost_(rhs.lastUpCost_),
180   lastDownDecrease_(rhs.lastDownDecrease_),
181   lastUpDecrease_(rhs.lastUpDecrease_),
182   numberTimesDown_(rhs.numberTimesDown_),
183   numberTimesUp_(rhs.numberTimesUp_),
184   numberTimesDownInfeasible_(rhs.numberTimesDownInfeasible_),
185   numberTimesUpInfeasible_(rhs.numberTimesUpInfeasible_),
186   numberBeforeTrust_(rhs.numberBeforeTrust_),
187   method_(rhs.method_)
188
189{
190}
191
192// Clone
193CbcObject *
194CbcSimpleIntegerDynamicPseudoCost::clone() const
195{
196  return new CbcSimpleIntegerDynamicPseudoCost(*this);
197}
198
199// Assignment operator
200CbcSimpleIntegerDynamicPseudoCost & 
201CbcSimpleIntegerDynamicPseudoCost::operator=( const CbcSimpleIntegerDynamicPseudoCost& rhs)
202{
203  if (this!=&rhs) {
204    CbcSimpleInteger::operator=(rhs);
205    downDynamicPseudoCost_=rhs.downDynamicPseudoCost_;
206    upDynamicPseudoCost_=rhs.upDynamicPseudoCost_;
207    upDownSeparator_=rhs.upDownSeparator_;
208    sumDownCost_ = rhs.sumDownCost_;
209    sumUpCost_ = rhs.sumUpCost_;
210    sumDownChange_ = rhs.sumDownChange_;
211    sumUpChange_ = rhs.sumUpChange_;
212    sumDownCostSquared_ = rhs.sumDownCostSquared_;
213    sumUpCostSquared_ = rhs.sumUpCostSquared_;
214    sumDownDecrease_ = rhs.sumDownDecrease_;
215    sumUpDecrease_ = rhs.sumUpDecrease_;
216    lastDownCost_ = rhs.lastDownCost_;
217    lastUpCost_ = rhs.lastUpCost_;
218    lastDownDecrease_ = rhs.lastDownDecrease_;
219    lastUpDecrease_ = rhs.lastUpDecrease_;
220    numberTimesDown_ = rhs.numberTimesDown_;
221    numberTimesUp_ = rhs.numberTimesUp_;
222    numberTimesDownInfeasible_ = rhs.numberTimesDownInfeasible_;
223    numberTimesUpInfeasible_ = rhs.numberTimesUpInfeasible_;
224    numberBeforeTrust_ = rhs.numberBeforeTrust_;
225    method_=rhs.method_;
226  }
227  return *this;
228}
229
230// Destructor
231CbcSimpleIntegerDynamicPseudoCost::~CbcSimpleIntegerDynamicPseudoCost ()
232{
233}
234// Creates a branching object
235CbcBranchingObject * 
236CbcSimpleIntegerDynamicPseudoCost::createBranch(int way) 
237{
238  const double * solution = model_->testSolution();
239  const double * lower = model_->getCbcColLower();
240  const double * upper = model_->getCbcColUpper();
241  double value = solution[columnNumber_];
242  value = CoinMax(value, lower[columnNumber_]);
243  value = CoinMin(value, upper[columnNumber_]);
244#ifndef NDEBUG
245  double nearest = floor(value+0.5);
246  double integerTolerance = 
247    model_->getDblParam(CbcModel::CbcIntegerTolerance);
248  assert (upper[columnNumber_]>lower[columnNumber_]);
249#endif
250  if (!model_->hotstartSolution()) {
251    assert (fabs(value-nearest)>integerTolerance);
252  } else {
253    const double * hotstartSolution = model_->hotstartSolution();
254    double targetValue = hotstartSolution[columnNumber_];
255    if (way>0)
256      value = targetValue-0.1;
257    else
258      value = targetValue+0.1;
259  }
260  CbcDynamicPseudoCostBranchingObject * newObject = 
261    new CbcDynamicPseudoCostBranchingObject(model_,sequence_,way,
262                                            value,this);
263  double up =  upDynamicPseudoCost_*(ceil(value)-value);
264  double down =  downDynamicPseudoCost_*(value-floor(value));
265  double changeInGuessed=up-down;
266  if (way>0)
267    changeInGuessed = - changeInGuessed;
268  changeInGuessed=CoinMax(0.0,changeInGuessed);
269  //if (way>0)
270  //changeInGuessed += 1.0e8; // bias to stay up
271  newObject->setChangeInGuessed(changeInGuessed);
272  newObject->setOriginalObject(this);
273  return newObject;
274}
275/* Create an OsiSolverBranch object
276   
277This returns NULL if branch not represented by bound changes
278*/
279OsiSolverBranch * 
280CbcSimpleIntegerDynamicPseudoCost::solverBranch() const
281{
282  OsiSolverInterface * solver = model_->solver();
283  const double * solution = model_->testSolution();
284  const double * lower = solver->getColLower();
285  const double * upper = solver->getColUpper();
286  double value = solution[columnNumber_];
287  value = CoinMax(value, lower[columnNumber_]);
288  value = CoinMin(value, upper[columnNumber_]);
289  assert (upper[columnNumber_]>lower[columnNumber_]);
290#ifndef NDEBUG
291  double nearest = floor(value+0.5);
292  double integerTolerance = 
293    model_->getDblParam(CbcModel::CbcIntegerTolerance);
294  assert (fabs(value-nearest)>integerTolerance);
295#endif
296  OsiSolverBranch * branch = new OsiSolverBranch();
297  branch->addBranch(columnNumber_,value);
298  return branch;
299}
300 
301// Infeasibility - large is 0.5
302double 
303CbcSimpleIntegerDynamicPseudoCost::infeasibility(int & preferredWay) const
304{
305  const double * solution = model_->testSolution();
306  const double * lower = model_->getCbcColLower();
307  const double * upper = model_->getCbcColUpper();
308  if (upper[columnNumber_]==lower[columnNumber_]) {
309    // fixed
310    preferredWay=1;
311    return 0.0;
312  }
313  assert (breakEven_>0.0&&breakEven_<1.0);
314  double value = solution[columnNumber_];
315  value = CoinMax(value, lower[columnNumber_]);
316  value = CoinMin(value, upper[columnNumber_]);
317  /*printf("%d %g %g %g %g\n",columnNumber_,value,lower[columnNumber_],
318    solution[columnNumber_],upper[columnNumber_]);*/
319  double nearest = floor(value+0.5);
320  double integerTolerance = 
321    model_->getDblParam(CbcModel::CbcIntegerTolerance);
322  double below = floor(value+integerTolerance);
323  double above = below+1.0;
324  if (above>upper[columnNumber_]) {
325    above=below;
326    below = above -1;
327  }
328#define INFEAS
329#ifdef INFEAS
330  double distanceToCutoff=0.0;
331  double objectiveValue = model_->getCurrentMinimizationObjValue();
332  distanceToCutoff =  model_->getCutoff()  - objectiveValue;
333  if (distanceToCutoff<1.0e20) 
334    distanceToCutoff *= 10.0;
335  else 
336    distanceToCutoff = 1.0e2 + fabs(objectiveValue);
337#endif
338  double sum;
339  int number;
340  double downCost = CoinMax(value-below,0.0);
341  sum = sumDownCost_;
342  number = numberTimesDown_;
343#ifdef INFEAS
344  sum += numberTimesDownInfeasible_*(distanceToCutoff/(downCost+1.0e-12));
345#endif
346  if (number>0)
347    downCost *= sum / (double) number;
348  else
349    downCost  *=  downDynamicPseudoCost_;
350  double upCost = CoinMax((above-value),0.0);
351  sum = sumUpCost_;
352  number = numberTimesUp_;
353#ifdef INFEAS
354  sum += numberTimesUpInfeasible_*(distanceToCutoff/(upCost+1.0e-12));
355#endif
356  if (number>0)
357    upCost *= sum / (double) number;
358  else
359    upCost  *=  upDynamicPseudoCost_;
360  if (downCost>=upCost)
361    preferredWay=1;
362  else
363    preferredWay=-1;
364  // See if up down choice set
365  if (upDownSeparator_>0.0) {
366    preferredWay = (value-below>=upDownSeparator_) ? 1 : -1;
367  }
368  if (preferredWay_)
369    preferredWay=preferredWay_;
370  // weight at 1.0 is max min
371#define WEIGHT_AFTER 0.9
372#define WEIGHT_BEFORE 0.3
373  if (fabs(value-nearest)<=integerTolerance) {
374    return 0.0;
375  } else {
376    int stateOfSearch = model_->stateOfSearch();
377    double returnValue=0.0;
378    double minValue = CoinMin(downCost,upCost);
379    double maxValue = CoinMax(downCost,upCost);
380    if (stateOfSearch<=1||model_->currentNode()->depth()<=10/* was ||maxValue>0.2*distanceToCutoff*/) {
381      // no solution
382      returnValue = WEIGHT_BEFORE*minValue + (1.0-WEIGHT_BEFORE)*maxValue;
383    } else {
384      // some solution
385      returnValue = WEIGHT_AFTER*minValue + (1.0-WEIGHT_AFTER)*maxValue;
386    }
387    if (numberTimesUp_<numberBeforeTrust_||
388        numberTimesDown_<numberBeforeTrust_) {
389      //if (returnValue<1.0e10)
390      //returnValue += 1.0e12;
391      //else
392      returnValue *= 1.0e3;
393      if (!numberTimesUp_&&!numberTimesDown_)
394        returnValue=1.0e50;
395    }
396    //if (fabs(value-0.5)<1.0e-5) {
397    //returnValue = 3.0*returnValue + 0.2;
398    //} else if (value>0.9) {
399    //returnValue = 2.0*returnValue + 0.1;
400    //}
401    return CoinMax(returnValue,1.0e-15);
402  }
403}
404
405// Return "up" estimate
406double 
407CbcSimpleIntegerDynamicPseudoCost::upEstimate() const
408{
409  const double * solution = model_->testSolution();
410  const double * lower = model_->getCbcColLower();
411  const double * upper = model_->getCbcColUpper();
412  double value = solution[columnNumber_];
413  value = CoinMax(value, lower[columnNumber_]);
414  value = CoinMin(value, upper[columnNumber_]);
415  if (upper[columnNumber_]==lower[columnNumber_]) {
416    // fixed
417    return 0.0;
418  }
419  double integerTolerance = 
420    model_->getDblParam(CbcModel::CbcIntegerTolerance);
421  double below = floor(value+integerTolerance);
422  double above = below+1.0;
423  if (above>upper[columnNumber_]) {
424    above=below;
425    below = above -1;
426  }
427  double upCost = CoinMax((above-value)*upDynamicPseudoCost_,0.0);
428  return upCost;
429}
430// Return "down" estimate
431double 
432CbcSimpleIntegerDynamicPseudoCost::downEstimate() const
433{
434  const double * solution = model_->testSolution();
435  const double * lower = model_->getCbcColLower();
436  const double * upper = model_->getCbcColUpper();
437  double value = solution[columnNumber_];
438  value = CoinMax(value, lower[columnNumber_]);
439  value = CoinMin(value, upper[columnNumber_]);
440  if (upper[columnNumber_]==lower[columnNumber_]) {
441    // fixed
442    return 0.0;
443  }
444  double integerTolerance = 
445    model_->getDblParam(CbcModel::CbcIntegerTolerance);
446  double below = floor(value+integerTolerance);
447  double above = below+1.0;
448  if (above>upper[columnNumber_]) {
449    above=below;
450    below = above -1;
451  }
452  double downCost = CoinMax((value-below)*downDynamicPseudoCost_,0.0);
453  return downCost;
454}
455// Print
456void 
457CbcSimpleIntegerDynamicPseudoCost::print(int type,double value) const
458{
459  if (!type) {
460    double meanDown =0.0;
461    double devDown =0.0;
462    if (numberTimesDown_) {
463      meanDown = sumDownCost_/(double) numberTimesDown_;
464      devDown = meanDown*meanDown + sumDownCostSquared_ - 
465        2.0*meanDown*sumDownCost_;
466      if (devDown>=0.0)
467        devDown = sqrt(devDown);
468    }
469    double meanUp =0.0;
470    double devUp =0.0;
471    if (numberTimesUp_) {
472      meanUp = sumUpCost_/(double) numberTimesUp_;
473      devUp = meanUp*meanUp + sumUpCostSquared_ - 
474        2.0*meanUp*sumUpCost_;
475      if (devUp>=0.0)
476        devUp = sqrt(devUp);
477    }
478#if 0
479    printf("%d down %d times (%d inf) mean %g (dev %g) up %d times (%d inf) mean %g (dev %g)\n",
480           columnNumber_,
481           numberTimesDown_,numberTimesDownInfeasible_,meanDown,devDown,
482           numberTimesUp_,numberTimesUpInfeasible_,meanUp,devUp);
483#else
484    printf("%d down %d times (%d inf) mean %g  up %d times (%d inf) mean %g\n",
485           columnNumber_,
486           numberTimesDown_,numberTimesDownInfeasible_,meanDown,
487           numberTimesUp_,numberTimesUpInfeasible_,meanUp);
488#endif
489  } else {
490    const double * upper = model_->getCbcColUpper();
491    double integerTolerance = 
492      model_->getDblParam(CbcModel::CbcIntegerTolerance);
493    double below = floor(value+integerTolerance);
494    double above = below+1.0;
495    if (above>upper[columnNumber_]) {
496      above=below;
497      below = above -1;
498    }
499    double objectiveValue = model_->getCurrentMinimizationObjValue();
500    double distanceToCutoff =  model_->getCutoff() - objectiveValue;
501    if (distanceToCutoff<1.0e20) 
502      distanceToCutoff *= 10.0;
503    else 
504      distanceToCutoff = 1.0e2 + fabs(objectiveValue);
505    double sum;
506    int number;
507    double downCost = CoinMax(value-below,0.0);
508    double downCost0 = downCost*downDynamicPseudoCost_;
509    sum = sumDownCost();
510    number = numberTimesDown();
511    sum += numberTimesDownInfeasible()*(distanceToCutoff/(downCost+1.0e-12));
512    if (number>0)
513      downCost *= sum / (double) number;
514    else
515      downCost  *=  downDynamicPseudoCost_;
516    double upCost = CoinMax((above-value),0.0);
517    double upCost0 = upCost*upDynamicPseudoCost_;
518    sum = sumUpCost();
519    number = numberTimesUp();
520    sum += numberTimesUpInfeasible()*(distanceToCutoff/(upCost+1.0e-12));
521    if (number>0)
522      upCost *= sum / (double) number;
523    else
524      upCost  *=  upDynamicPseudoCost_;
525    printf("%d down %d times %g (est %g)  up %d times %g (est %g)\n",
526           columnNumber_,
527           numberTimesDown_,downCost,downCost0,
528           numberTimesUp_,upCost,upCost0);
529  }
530}
531
532// Default Constructor
533CbcDynamicPseudoCostBranchingObject::CbcDynamicPseudoCostBranchingObject()
534  :CbcIntegerBranchingObject()
535{
536  changeInGuessed_=1.0e-5;
537  object_=NULL;
538}
539
540// Useful constructor
541CbcDynamicPseudoCostBranchingObject::CbcDynamicPseudoCostBranchingObject (CbcModel * model, 
542                                                                          int variable, 
543                                                                          int way , double value,
544                                       CbcSimpleIntegerDynamicPseudoCost * object) 
545  :CbcIntegerBranchingObject(model,variable,way,value)
546{
547  changeInGuessed_=1.0e-5;
548  object_=object;
549}
550// Useful constructor for fixing
551CbcDynamicPseudoCostBranchingObject::CbcDynamicPseudoCostBranchingObject (CbcModel * model, 
552                                                      int variable, int way,
553                                                      double lowerValue, 
554                                                      double upperValue)
555  :CbcIntegerBranchingObject(model,variable,way,lowerValue)
556{
557  changeInGuessed_=1.0e100;
558  object_=NULL;
559}
560 
561
562// Copy constructor
563CbcDynamicPseudoCostBranchingObject::CbcDynamicPseudoCostBranchingObject ( 
564                                 const CbcDynamicPseudoCostBranchingObject & rhs)
565  :CbcIntegerBranchingObject(rhs)
566{
567  changeInGuessed_ = rhs.changeInGuessed_;
568  object_=rhs.object_;
569}
570
571// Assignment operator
572CbcDynamicPseudoCostBranchingObject & 
573CbcDynamicPseudoCostBranchingObject::operator=( const CbcDynamicPseudoCostBranchingObject& rhs)
574{
575  if (this != &rhs) {
576    CbcIntegerBranchingObject::operator=(rhs);
577    changeInGuessed_ = rhs.changeInGuessed_;
578    object_=rhs.object_;
579  }
580  return *this;
581}
582CbcBranchingObject * 
583CbcDynamicPseudoCostBranchingObject::clone() const
584{ 
585  return (new CbcDynamicPseudoCostBranchingObject(*this));
586}
587
588
589// Destructor
590CbcDynamicPseudoCostBranchingObject::~CbcDynamicPseudoCostBranchingObject ()
591{
592}
593
594/*
595  Perform a branch by adjusting the bounds of the specified variable. Note
596  that each arm of the branch advances the object to the next arm by
597  advancing the value of way_.
598
599  Providing new values for the variable's lower and upper bounds for each
600  branching direction gives a little bit of additional flexibility and will
601  be easily extensible to multi-way branching.
602  Returns change in guessed objective on next branch
603*/
604double
605CbcDynamicPseudoCostBranchingObject::branch(bool normalBranch)
606{
607  CbcIntegerBranchingObject::branch(normalBranch);
608  return changeInGuessed_;
609}
610/* Some branchingObjects may claim to be able to skip
611   strong branching.  If so they have to fill in CbcStrongInfo.
612   The object mention in incoming CbcStrongInfo must match.
613   Returns nonzero if skip is wanted */
614int 
615CbcDynamicPseudoCostBranchingObject::fillStrongInfo( CbcStrongInfo & info)
616{
617  assert (object_);
618  assert (info.possibleBranch==this);
619  if (object_->numberTimesUp()<object_->numberBeforeTrust()||
620      object_->numberTimesDown()<object_->numberBeforeTrust()) {
621    return 0;
622  } else {
623    info.upMovement = object_->upDynamicPseudoCost()*(ceil(value_)-value_);
624    info.downMovement = object_->downDynamicPseudoCost()*(value_-floor(value_));
625    info.numIntInfeasUp  -= (int) (object_->sumUpDecrease()/
626                                   ((double) object_->numberTimesUp()));
627    info.numObjInfeasUp = 0;
628    info.finishedUp = false;
629    info.numItersUp = 0;
630    info.numIntInfeasDown  -= (int) (object_->sumDownDecrease()/
631                                   ((double) object_->numberTimesDown()));
632    info.numObjInfeasDown = 0;
633    info.finishedDown = false;
634    info.numItersDown = 0;
635    info.fix =0;
636    return 1;
637  }
638}
639 
640// Default Constructor
641CbcBranchDynamicDecision::CbcBranchDynamicDecision()
642  :CbcBranchDecision()
643{
644  bestCriterion_ = 0.0;
645  bestChangeUp_ = 0.0;
646  bestNumberUp_ = 0;
647  bestChangeDown_ = 0.0;
648  bestNumberDown_ = 0;
649  bestObject_ = NULL;
650}
651
652// Copy constructor
653CbcBranchDynamicDecision::CbcBranchDynamicDecision (
654                                    const CbcBranchDynamicDecision & rhs)
655  :CbcBranchDecision()
656{
657  bestCriterion_ = rhs.bestCriterion_;
658  bestChangeUp_ = rhs.bestChangeUp_;
659  bestNumberUp_ = rhs.bestNumberUp_;
660  bestChangeDown_ = rhs.bestChangeDown_;
661  bestNumberDown_ = rhs.bestNumberDown_;
662  bestObject_ = rhs.bestObject_;
663}
664
665CbcBranchDynamicDecision::~CbcBranchDynamicDecision()
666{
667}
668
669// Clone
670CbcBranchDecision * 
671CbcBranchDynamicDecision::clone() const
672{
673  return new CbcBranchDynamicDecision(*this);
674}
675
676// Initialize i.e. before start of choosing at a node
677void 
678CbcBranchDynamicDecision::initialize(CbcModel * model)
679{
680  bestCriterion_ = 0.0;
681  bestChangeUp_ = 0.0;
682  bestNumberUp_ = 0;
683  bestChangeDown_ = 0.0;
684  bestNumberDown_ = 0;
685  bestObject_ = NULL;
686}
687
688/* Saves a clone of current branching object.  Can be used to update
689      information on object causing branch - after branch */
690void 
691CbcBranchDynamicDecision::saveBranchingObject(CbcBranchingObject * object) 
692{
693  object_ = object->clone();
694}
695/* Pass in information on branch just done.
696   assumes object can get information from solver */
697void 
698CbcBranchDynamicDecision::updateInformation(OsiSolverInterface * solver,
699                                            const CbcNode * node)
700{
701  assert (object_);
702  const CbcModel * model = object_->model();
703  double originalValue=node->objectiveValue();
704  int originalUnsatisfied = node->numberUnsatisfied();
705  double objectiveValue = solver->getObjValue()*model->getObjSense();
706  int unsatisfied=0;
707  int i;
708  int numberIntegers = model->numberIntegers();;
709  const double * solution = solver->getColSolution();
710  //const double * lower = solver->getColLower();
711  //const double * upper = solver->getColUpper();
712  CbcDynamicPseudoCostBranchingObject * branchingObject =
713    dynamic_cast<CbcDynamicPseudoCostBranchingObject *>(object_);
714  assert (branchingObject);
715  CbcSimpleIntegerDynamicPseudoCost *  object = branchingObject->object();
716  double change = CoinMax(0.0,objectiveValue-originalValue);
717  // probably should also ignore if stopped
718  int iStatus;
719  if (solver->isProvenOptimal())
720    iStatus=0; // optimal
721  else if (solver->isIterationLimitReached()
722           &&!solver->isDualObjectiveLimitReached())
723    iStatus=2; // unknown
724  else
725    iStatus=1; // infeasible
726
727  bool feasible = iStatus!=1;
728  if (feasible) {
729    double integerTolerance = 
730      model->getDblParam(CbcModel::CbcIntegerTolerance);
731    const int * integerVariable = model->integerVariable();
732    for (i=0;i<numberIntegers;i++) {
733      int j=integerVariable[i];
734      double value = solution[j];
735      double nearest = floor(value+0.5);
736      if (fabs(value-nearest)>integerTolerance) 
737        unsatisfied++;
738    }
739  }
740  int way = object_->way();
741  double value = object_->value();
742#define TYPE2
743  if (way<0) {
744    // down
745    if (feasible) {
746      object->incrementNumberTimesDown();
747      object->addToSumDownChange(1.0e-30+value-floor(value));
748      object->addToSumDownDecrease(originalUnsatisfied-unsatisfied);
749#ifndef TYPE2
750      object->addToSumDownCost(change/(1.0e-30+(value-floor(value))));
751      object->setDownDynamicPseudoCost(object->sumDownCost()/(double) object->numberTimesDown());
752#else
753      object->addToSumDownCost(change);
754      object->setDownDynamicPseudoCost(object->sumDownCost()/object->sumDownChange());
755#endif
756    } else {
757      object->incrementNumberTimesDownInfeasible();
758    }
759  } else {
760    // up
761    if (feasible) {
762      object->incrementNumberTimesUp();
763      object->addToSumUpChange(1.0e-30+ceil(value)-value);
764      object->addToSumUpDecrease(unsatisfied-originalUnsatisfied);
765#ifndef TYPE2
766      object->addToSumUpCost(change/(1.0e-30+(ceil(value)-value)));
767      object->setUpDynamicPseudoCost(object->sumUpCost()/(double) object->numberTimesUp());
768#else
769      object->addToSumUpCost(change);
770      object->setUpDynamicPseudoCost(object->sumUpCost()/object->sumUpChange());
771#endif
772    } else {
773      object->incrementNumberTimesUpInfeasible();
774    }
775  }
776  delete object_;
777  object_=NULL;
778}
779
780/*
781  Simple dynamic decision algorithm. Compare based on infeasibility (numInfUp,
782  numInfDown) until a solution is found by search, then switch to change in
783  objective (changeUp, changeDown). Note that bestSoFar is remembered in
784  bestObject_, so the parameter bestSoFar is unused.
785*/
786
787int
788CbcBranchDynamicDecision::betterBranch(CbcBranchingObject * thisOne,
789                            CbcBranchingObject * bestSoFar,
790                            double changeUp, int numInfUp,
791                            double changeDown, int numInfDown)
792{
793  int stateOfSearch = thisOne->model()->stateOfSearch();
794  int betterWay=0;
795  double value=0.0;
796  if (stateOfSearch<=1&&thisOne->model()->currentNode()->depth()>10) {
797#if 0
798    if (!bestObject_) {
799      bestNumberUp_=INT_MAX;
800      bestNumberDown_=INT_MAX;
801    }
802    // before solution - choose smallest number
803    // could add in depth as well
804    int bestNumber = CoinMin(bestNumberUp_,bestNumberDown_);
805    if (numInfUp<numInfDown) {
806      if (numInfUp<bestNumber) {
807        betterWay = 1;
808      } else if (numInfUp==bestNumber) {
809        if (changeUp<bestChangeUp_)
810          betterWay=1;
811      }
812    } else if (numInfUp>numInfDown) {
813      if (numInfDown<bestNumber) {
814        betterWay = -1;
815      } else if (numInfDown==bestNumber) {
816        if (changeDown<bestChangeDown_)
817          betterWay=-1;
818      }
819    } else {
820      // up and down have same number
821      bool better=false;
822      if (numInfUp<bestNumber) {
823        better=true;
824      } else if (numInfUp==bestNumber) {
825        if (CoinMin(changeUp,changeDown)<CoinMin(bestChangeUp_,bestChangeDown_)-1.0e-5)
826          better=true;;
827      }
828      if (better) {
829        // see which way
830        if (changeUp<=changeDown)
831          betterWay=1;
832        else
833          betterWay=-1;
834      }
835    }
836    if (betterWay) {
837      value = CoinMin(numInfUp,numInfDown);
838    }
839#else
840    if (!bestObject_) {
841      bestCriterion_=-1.0;
842    }
843    // got a solution
844    double minValue = CoinMin(changeDown,changeUp);
845    double maxValue = CoinMax(changeDown,changeUp);
846    value = WEIGHT_BEFORE*minValue + (1.0-WEIGHT_BEFORE)*maxValue;
847    if (value>bestCriterion_+1.0e-8) {
848      if (changeUp<=changeDown) {
849        betterWay=1;
850      } else {
851        betterWay=-1;
852      }
853    }
854#endif
855  } else {
856    if (!bestObject_) {
857      bestCriterion_=-1.0;
858    }
859    // got a solution
860    double minValue = CoinMin(changeDown,changeUp);
861    double maxValue = CoinMax(changeDown,changeUp);
862    // Reduce
863    maxValue = CoinMin(maxValue,minValue*2.0);
864    value = WEIGHT_AFTER*minValue + (1.0-WEIGHT_AFTER)*maxValue;
865    if (value>bestCriterion_+1.0e-8) {
866      if (changeUp<=changeDown) {
867        betterWay=1;
868      } else {
869        betterWay=-1;
870      }
871    }
872  }
873  if (betterWay) {
874    // maybe change better way
875    CbcDynamicPseudoCostBranchingObject * branchingObject =
876      dynamic_cast<CbcDynamicPseudoCostBranchingObject *>(thisOne);
877    assert (branchingObject);
878    CbcSimpleIntegerDynamicPseudoCost *  object = branchingObject->object();
879    double separator = object->upDownSeparator();
880    if (separator>0.0) {
881      const double * solution = thisOne->model()->testSolution();
882      double valueVariable = solution[object->columnNumber()];
883      betterWay = (valueVariable-floor(valueVariable)>=separator) ? 1 : -1;
884    }
885    bestCriterion_ = value;
886    bestChangeUp_ = changeUp;
887    bestNumberUp_ = numInfUp;
888    bestChangeDown_ = changeDown;
889    bestNumberDown_ = numInfDown;
890    bestObject_=thisOne;
891    // See if user is overriding way
892    if (thisOne->object()&&thisOne->object()->preferredWay())
893      betterWay = thisOne->object()->preferredWay();
894  }
895  return betterWay;
896}
897/* Sets or gets best criterion so far */
898void 
899CbcBranchDynamicDecision::setBestCriterion(double value)
900{ 
901  bestCriterion_ = value;
902}
903double 
904CbcBranchDynamicDecision::getBestCriterion() const
905{ 
906  return bestCriterion_;
907}
Note: See TracBrowser for help on using the repository browser.