source: trunk/CbcBranchDynamic.cpp @ 147

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

dynamic branching stuff

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