source: branches/devel/Cbc/src/CbcBranchDynamic.cpp @ 435

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

many changes

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