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

Last change on this file since 662 was 662, checked in by forrest, 12 years ago

faster parallel?

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 46.0 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,
60                                    int iColumn, double breakEven)
61  : CbcSimpleInteger(model,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,
123                                    int iColumn, double downDynamicPseudoCost,
124                                                        double upDynamicPseudoCost)
125  : CbcSimpleInteger(model,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/** Useful constructor
179
180  Loads dynamic upper & lower bounds for the specified variable.
181*/
182CbcSimpleIntegerDynamicPseudoCost::CbcSimpleIntegerDynamicPseudoCost (CbcModel * model,
183                                    int dummy, int iColumn, double downDynamicPseudoCost,
184                                                        double upDynamicPseudoCost)
185{
186  CbcSimpleIntegerDynamicPseudoCost(model,iColumn,downDynamicPseudoCost,upDynamicPseudoCost);
187}
188
189// Copy constructor
190CbcSimpleIntegerDynamicPseudoCost::CbcSimpleIntegerDynamicPseudoCost ( const CbcSimpleIntegerDynamicPseudoCost & rhs)
191  :CbcSimpleInteger(rhs),
192   downDynamicPseudoCost_(rhs.downDynamicPseudoCost_),
193   upDynamicPseudoCost_(rhs.upDynamicPseudoCost_),
194   upDownSeparator_(rhs.upDownSeparator_),
195   sumDownCost_(rhs.sumDownCost_),
196   sumUpCost_(rhs.sumUpCost_),
197   sumDownChange_(rhs.sumDownChange_),
198   sumUpChange_(rhs.sumUpChange_),
199   sumDownCostSquared_(rhs.sumDownCostSquared_),
200   sumUpCostSquared_(rhs.sumUpCostSquared_),
201   sumDownDecrease_(rhs.sumDownDecrease_),
202   sumUpDecrease_(rhs.sumUpDecrease_),
203   lastDownCost_(rhs.lastDownCost_),
204   lastUpCost_(rhs.lastUpCost_),
205   lastDownDecrease_(rhs.lastDownDecrease_),
206   lastUpDecrease_(rhs.lastUpDecrease_),
207   numberTimesDown_(rhs.numberTimesDown_),
208   numberTimesUp_(rhs.numberTimesUp_),
209   numberTimesDownInfeasible_(rhs.numberTimesDownInfeasible_),
210   numberTimesUpInfeasible_(rhs.numberTimesUpInfeasible_),
211   numberBeforeTrust_(rhs.numberBeforeTrust_),
212   numberTimesDownLocalFixed_(rhs.numberTimesDownLocalFixed_),
213   numberTimesUpLocalFixed_(rhs.numberTimesUpLocalFixed_),
214   numberTimesDownTotalFixed_(rhs.numberTimesDownTotalFixed_),
215   numberTimesUpTotalFixed_(rhs.numberTimesUpTotalFixed_),
216   numberTimesProbingTotal_(rhs.numberTimesProbingTotal_),
217   method_(rhs.method_)
218
219{
220}
221
222// Clone
223CbcObject *
224CbcSimpleIntegerDynamicPseudoCost::clone() const
225{
226  return new CbcSimpleIntegerDynamicPseudoCost(*this);
227}
228
229// Assignment operator
230CbcSimpleIntegerDynamicPseudoCost & 
231CbcSimpleIntegerDynamicPseudoCost::operator=( const CbcSimpleIntegerDynamicPseudoCost& rhs)
232{
233  if (this!=&rhs) {
234    CbcSimpleInteger::operator=(rhs);
235    downDynamicPseudoCost_=rhs.downDynamicPseudoCost_;
236    upDynamicPseudoCost_=rhs.upDynamicPseudoCost_;
237    upDownSeparator_=rhs.upDownSeparator_;
238    sumDownCost_ = rhs.sumDownCost_;
239    sumUpCost_ = rhs.sumUpCost_;
240    sumDownChange_ = rhs.sumDownChange_;
241    sumUpChange_ = rhs.sumUpChange_;
242    sumDownCostSquared_ = rhs.sumDownCostSquared_;
243    sumUpCostSquared_ = rhs.sumUpCostSquared_;
244    sumDownDecrease_ = rhs.sumDownDecrease_;
245    sumUpDecrease_ = rhs.sumUpDecrease_;
246    lastDownCost_ = rhs.lastDownCost_;
247    lastUpCost_ = rhs.lastUpCost_;
248    lastDownDecrease_ = rhs.lastDownDecrease_;
249    lastUpDecrease_ = rhs.lastUpDecrease_;
250    numberTimesDown_ = rhs.numberTimesDown_;
251    numberTimesUp_ = rhs.numberTimesUp_;
252    numberTimesDownInfeasible_ = rhs.numberTimesDownInfeasible_;
253    numberTimesUpInfeasible_ = rhs.numberTimesUpInfeasible_;
254    numberBeforeTrust_ = rhs.numberBeforeTrust_;
255    numberTimesDownLocalFixed_ = rhs.numberTimesDownLocalFixed_;
256    numberTimesUpLocalFixed_ = rhs.numberTimesUpLocalFixed_;
257    numberTimesDownTotalFixed_ = rhs.numberTimesDownTotalFixed_;
258    numberTimesUpTotalFixed_ = rhs.numberTimesUpTotalFixed_;
259    numberTimesProbingTotal_ = rhs.numberTimesProbingTotal_;
260    method_=rhs.method_;
261  }
262  return *this;
263}
264
265// Destructor
266CbcSimpleIntegerDynamicPseudoCost::~CbcSimpleIntegerDynamicPseudoCost ()
267{
268}
269// Copy some information i.e. just variable stuff
270void 
271CbcSimpleIntegerDynamicPseudoCost::copySome(CbcSimpleIntegerDynamicPseudoCost * otherObject)
272{
273  downDynamicPseudoCost_=otherObject->downDynamicPseudoCost_;
274  upDynamicPseudoCost_=otherObject->upDynamicPseudoCost_;
275  sumDownCost_ = otherObject->sumDownCost_;
276  sumUpCost_ = otherObject->sumUpCost_;
277  sumDownChange_ = otherObject->sumDownChange_;
278  sumUpChange_ = otherObject->sumUpChange_;
279  sumDownCostSquared_ = otherObject->sumDownCostSquared_;
280  sumUpCostSquared_ = otherObject->sumUpCostSquared_;
281  sumDownDecrease_ = otherObject->sumDownDecrease_;
282  sumUpDecrease_ = otherObject->sumUpDecrease_;
283  lastDownCost_ = otherObject->lastDownCost_;
284  lastUpCost_ = otherObject->lastUpCost_;
285  lastDownDecrease_ = otherObject->lastDownDecrease_;
286  lastUpDecrease_ = otherObject->lastUpDecrease_;
287  numberTimesDown_ = otherObject->numberTimesDown_;
288  numberTimesUp_ = otherObject->numberTimesUp_;
289  numberTimesDownInfeasible_ = otherObject->numberTimesDownInfeasible_;
290  numberTimesUpInfeasible_ = otherObject->numberTimesUpInfeasible_;
291  numberTimesDownLocalFixed_ = otherObject->numberTimesDownLocalFixed_;
292  numberTimesUpLocalFixed_ = otherObject->numberTimesUpLocalFixed_;
293  numberTimesDownTotalFixed_ = otherObject->numberTimesDownTotalFixed_;
294  numberTimesUpTotalFixed_ = otherObject->numberTimesUpTotalFixed_;
295  numberTimesProbingTotal_ = otherObject->numberTimesProbingTotal_;
296}
297// Creates a branching object
298CbcBranchingObject * 
299CbcSimpleIntegerDynamicPseudoCost::createBranch(int way) 
300{
301  const double * solution = model_->testSolution();
302  const double * lower = model_->getCbcColLower();
303  const double * upper = model_->getCbcColUpper();
304  double value = solution[columnNumber_];
305  value = CoinMax(value, lower[columnNumber_]);
306  value = CoinMin(value, upper[columnNumber_]);
307#ifndef NDEBUG
308  double nearest = floor(value+0.5);
309  double integerTolerance = 
310    model_->getDblParam(CbcModel::CbcIntegerTolerance);
311  assert (upper[columnNumber_]>lower[columnNumber_]);
312#endif
313  if (!model_->hotstartSolution()) {
314    assert (fabs(value-nearest)>integerTolerance);
315  } else {
316    const double * hotstartSolution = model_->hotstartSolution();
317    double targetValue = hotstartSolution[columnNumber_];
318    if (way>0)
319      value = targetValue-0.1;
320    else
321      value = targetValue+0.1;
322  }
323  CbcDynamicPseudoCostBranchingObject * newObject = 
324    new CbcDynamicPseudoCostBranchingObject(model_,columnNumber_,way,
325                                            value,this);
326  double up =  upDynamicPseudoCost_*(ceil(value)-value);
327  double down =  downDynamicPseudoCost_*(value-floor(value));
328  double changeInGuessed=up-down;
329  if (way>0)
330    changeInGuessed = - changeInGuessed;
331  changeInGuessed=CoinMax(0.0,changeInGuessed);
332  //if (way>0)
333  //changeInGuessed += 1.0e8; // bias to stay up
334  newObject->setChangeInGuessed(changeInGuessed);
335  newObject->setOriginalObject(this);
336  return newObject;
337}
338/* Create an OsiSolverBranch object
339   
340This returns NULL if branch not represented by bound changes
341*/
342OsiSolverBranch * 
343CbcSimpleIntegerDynamicPseudoCost::solverBranch() const
344{
345  OsiSolverInterface * solver = model_->solver();
346  const double * solution = model_->testSolution();
347  const double * lower = solver->getColLower();
348  const double * upper = solver->getColUpper();
349  double value = solution[columnNumber_];
350  value = CoinMax(value, lower[columnNumber_]);
351  value = CoinMin(value, upper[columnNumber_]);
352  assert (upper[columnNumber_]>lower[columnNumber_]);
353#ifndef NDEBUG
354  double nearest = floor(value+0.5);
355  double integerTolerance = 
356    model_->getDblParam(CbcModel::CbcIntegerTolerance);
357  assert (fabs(value-nearest)>integerTolerance);
358#endif
359  OsiSolverBranch * branch = new OsiSolverBranch();
360  branch->addBranch(columnNumber_,value);
361  return branch;
362}
363 
364// Infeasibility - large is 0.5
365double 
366CbcSimpleIntegerDynamicPseudoCost::infeasibility(int & preferredWay) const
367{
368  const double * solution = model_->testSolution();
369  const double * lower = model_->getCbcColLower();
370  const double * upper = model_->getCbcColUpper();
371  if (upper[columnNumber_]==lower[columnNumber_]) {
372    // fixed
373    preferredWay=1;
374    return 0.0;
375  }
376  assert (breakEven_>0.0&&breakEven_<1.0);
377  double value = solution[columnNumber_];
378  value = CoinMax(value, lower[columnNumber_]);
379  value = CoinMin(value, upper[columnNumber_]);
380  /*printf("%d %g %g %g %g\n",columnNumber_,value,lower[columnNumber_],
381    solution[columnNumber_],upper[columnNumber_]);*/
382  double nearest = floor(value+0.5);
383  double integerTolerance = 
384    model_->getDblParam(CbcModel::CbcIntegerTolerance);
385  double below = floor(value+integerTolerance);
386  double above = below+1.0;
387  if (above>upper[columnNumber_]) {
388    above=below;
389    below = above -1;
390  }
391#define INFEAS 1
392#if INFEAS==1
393  double distanceToCutoff=0.0;
394  double objectiveValue = model_->getCurrentMinimizationObjValue();
395  distanceToCutoff =  model_->getCutoff()  - objectiveValue;
396  if (distanceToCutoff<1.0e20) 
397    distanceToCutoff *= 10.0;
398  else 
399    distanceToCutoff = 1.0e2 + fabs(objectiveValue);
400#endif
401  double sum;
402  int number;
403  double downCost = CoinMax(value-below,0.0);
404  sum = sumDownCost_;
405  number = numberTimesDown_;
406#if INFEAS==1
407  sum += numberTimesDownInfeasible_*(distanceToCutoff/(downCost+1.0e-12));
408  number += numberTimesDownInfeasible_;
409#endif
410  if (number>0)
411    downCost *= sum / (double) number;
412  else
413    downCost  *=  downDynamicPseudoCost_;
414  double upCost = CoinMax((above-value),0.0);
415  sum = sumUpCost_;
416  number = numberTimesUp_;
417#if INFEAS==1
418  sum += numberTimesUpInfeasible_*(distanceToCutoff/(upCost+1.0e-12));
419  number += numberTimesUpInfeasible_;
420#endif
421  if (number>0)
422    upCost *= sum / (double) number;
423  else
424    upCost  *=  upDynamicPseudoCost_;
425  if (downCost>=upCost)
426    preferredWay=1;
427  else
428    preferredWay=-1;
429  // See if up down choice set
430  if (upDownSeparator_>0.0) {
431    preferredWay = (value-below>=upDownSeparator_) ? 1 : -1;
432  }
433  if (preferredWay_)
434    preferredWay=preferredWay_;
435  // weight at 1.0 is max min
436#define WEIGHT_AFTER 0.9
437#define WEIGHT_BEFORE 0.3
438  if (fabs(value-nearest)<=integerTolerance) {
439    return 0.0;
440  } else {
441    int stateOfSearch = model_->stateOfSearch();
442    double returnValue=0.0;
443    double minValue = CoinMin(downCost,upCost);
444    double maxValue = CoinMax(downCost,upCost);
445    if (stateOfSearch<=1||model_->currentNode()->depth()<=10/* was ||maxValue>0.2*distanceToCutoff*/) {
446      // no solution
447      returnValue = WEIGHT_BEFORE*minValue + (1.0-WEIGHT_BEFORE)*maxValue;
448    } else {
449      // some solution
450      returnValue = WEIGHT_AFTER*minValue + (1.0-WEIGHT_AFTER)*maxValue;
451    }
452    if (numberTimesUp_<numberBeforeTrust_||
453        numberTimesDown_<numberBeforeTrust_) {
454      //if (returnValue<1.0e10)
455      //returnValue += 1.0e12;
456      //else
457      returnValue *= 1.0e3;
458      if (!numberTimesUp_&&!numberTimesDown_)
459        returnValue=1.0e50;
460    }
461    //if (fabs(value-0.5)<1.0e-5) {
462    //returnValue = 3.0*returnValue + 0.2;
463    //} else if (value>0.9) {
464    //returnValue = 2.0*returnValue + 0.1;
465    //}
466    if (method_==1) {
467      // probing
468      // average
469      double up=1.0e-15;
470      double down=1.0e-15;
471      if (numberTimesProbingTotal_) {
472        up += numberTimesUpTotalFixed_/((double) numberTimesProbingTotal_);
473        down += numberTimesDownTotalFixed_/((double) numberTimesProbingTotal_);
474      }
475      returnValue = 1 + 10.0*CoinMin(numberTimesDownLocalFixed_,numberTimesUpLocalFixed_) +
476        CoinMin(down,up);
477      returnValue *= 1.0e-3;
478    }
479    return CoinMax(returnValue,1.0e-15);
480  }
481}
482
483double
484CbcSimpleIntegerDynamicPseudoCost::infeasibility(const OsiSolverInterface * solver, const OsiBranchingInformation * info,
485                         int & preferredWay) const
486{
487  double value = info->solution_[columnNumber_];
488  value = CoinMax(value, info->lower_[columnNumber_]);
489  value = CoinMin(value, info->upper_[columnNumber_]);
490  if (info->upper_[columnNumber_]==info->lower_[columnNumber_]) {
491    // fixed
492    preferredWay=1;
493    return 0.0;
494  }
495  assert (breakEven_>0.0&&breakEven_<1.0);
496  double nearest = floor(value+0.5);
497  double integerTolerance = info->integerTolerance_; 
498  double below = floor(value+integerTolerance);
499  double above = below+1.0;
500  if (above>info->upper_[columnNumber_]) {
501    above=below;
502    below = above -1;
503  }
504#if INFEAS==1
505  double objectiveValue = info->objectiveValue_;
506  double distanceToCutoff =  info->cutoff_  - objectiveValue;
507  if (distanceToCutoff<1.0e20) 
508    distanceToCutoff *= 10.0;
509  else 
510    distanceToCutoff = 1.0e2 + fabs(objectiveValue);
511#endif
512  double sum;
513  int number;
514  double downCost = CoinMax(value-below,0.0);
515  sum = sumDownCost_;
516  number = numberTimesDown_;
517#if INFEAS==1
518  sum += numberTimesDownInfeasible_*(distanceToCutoff/(downCost+1.0e-12));
519  number += numberTimesDownInfeasible_;
520#endif
521  if (number>0)
522    downCost *= sum / (double) number;
523  else
524    downCost  *=  downDynamicPseudoCost_;
525  double upCost = CoinMax((above-value),0.0);
526  sum = sumUpCost_;
527  number = numberTimesUp_;
528#if INFEAS==1
529  sum += numberTimesUpInfeasible_*(distanceToCutoff/(upCost+1.0e-12));
530  number += numberTimesUpInfeasible_;
531#endif
532  if (number>0)
533    upCost *= sum / (double) number;
534  else
535    upCost  *=  upDynamicPseudoCost_;
536  if (downCost>=upCost)
537    preferredWay=1;
538  else
539    preferredWay=-1;
540  // See if up down choice set
541  if (upDownSeparator_>0.0) {
542    preferredWay = (value-below>=upDownSeparator_) ? 1 : -1;
543  }
544  if (preferredWay_)
545    preferredWay=preferredWay_;
546  // weight at 1.0 is max min
547#define WEIGHT_BEFORE 0.3
548  if (fabs(value-nearest)<=integerTolerance) {
549    return 0.0;
550  } else {
551    double returnValue=0.0;
552    double minValue = CoinMin(downCost,upCost);
553    double maxValue = CoinMax(downCost,upCost);
554    if (!info->numberBranchingSolutions_||info->depth_<=10/* was ||maxValue>0.2*distanceToCutoff*/) {
555      // no solution
556      returnValue = WEIGHT_BEFORE*minValue + (1.0-WEIGHT_BEFORE)*maxValue;
557    } else {
558      // some solution
559      returnValue = WEIGHT_AFTER*minValue + (1.0-WEIGHT_AFTER)*maxValue;
560    }
561    if (numberTimesUp_<numberBeforeTrust_||
562        numberTimesDown_<numberBeforeTrust_) {
563      //if (returnValue<1.0e10)
564      //returnValue += 1.0e12;
565      //else
566      returnValue *= 1.0e3;
567      if (!numberTimesUp_&&!numberTimesDown_)
568        returnValue=1.0e50;
569    }
570    //if (fabs(value-0.5)<1.0e-5) {
571    //returnValue = 3.0*returnValue + 0.2;
572    //} else if (value>0.9) {
573    //returnValue = 2.0*returnValue + 0.1;
574    //}
575    if (method_==1) {
576      // probing
577      // average
578      double up=1.0e-15;
579      double down=1.0e-15;
580      if (numberTimesProbingTotal_) {
581        up += numberTimesUpTotalFixed_/((double) numberTimesProbingTotal_);
582        down += numberTimesDownTotalFixed_/((double) numberTimesProbingTotal_);
583      }
584      returnValue = 1 + 10.0*CoinMin(numberTimesDownLocalFixed_,numberTimesUpLocalFixed_) +
585        CoinMin(down,up);
586      returnValue *= 1.0e-3;
587    }
588    return CoinMax(returnValue,1.0e-15);
589  }
590}
591// Creates a branching object
592CbcBranchingObject * 
593CbcSimpleIntegerDynamicPseudoCost::createBranch(OsiSolverInterface * solver, const OsiBranchingInformation * info, int way) 
594{
595  double value = info->solution_[columnNumber_];
596  value = CoinMax(value, info->lower_[columnNumber_]);
597  value = CoinMin(value, info->upper_[columnNumber_]);
598  assert (info->upper_[columnNumber_]>info->lower_[columnNumber_]);
599  if (!info->hotstartSolution_) {
600#ifndef NDEBUG
601    double nearest = floor(value+0.5);
602    assert (fabs(value-nearest)>info->integerTolerance_);
603#endif
604  } else {
605    double targetValue = info->hotstartSolution_[columnNumber_];
606    if (way>0)
607      value = targetValue-0.1;
608    else
609      value = targetValue+0.1;
610  }
611  CbcDynamicPseudoCostBranchingObject * newObject = 
612    new CbcDynamicPseudoCostBranchingObject(model_,columnNumber_,way,
613                                            value,this);
614  double up =  upDynamicPseudoCost_*(ceil(value)-value);
615  double down =  downDynamicPseudoCost_*(value-floor(value));
616  double changeInGuessed=up-down;
617  if (way>0)
618    changeInGuessed = - changeInGuessed;
619  changeInGuessed=CoinMax(0.0,changeInGuessed);
620  //if (way>0)
621  //changeInGuessed += 1.0e8; // bias to stay up
622  newObject->setChangeInGuessed(changeInGuessed);
623  newObject->setOriginalObject(this);
624  return newObject;
625}
626
627// Return "up" estimate
628double 
629CbcSimpleIntegerDynamicPseudoCost::upEstimate() const
630{
631  const double * solution = model_->testSolution();
632  const double * lower = model_->getCbcColLower();
633  const double * upper = model_->getCbcColUpper();
634  double value = solution[columnNumber_];
635  value = CoinMax(value, lower[columnNumber_]);
636  value = CoinMin(value, upper[columnNumber_]);
637  if (upper[columnNumber_]==lower[columnNumber_]) {
638    // fixed
639    return 0.0;
640  }
641  double integerTolerance = 
642    model_->getDblParam(CbcModel::CbcIntegerTolerance);
643  double below = floor(value+integerTolerance);
644  double above = below+1.0;
645  if (above>upper[columnNumber_]) {
646    above=below;
647    below = above -1;
648  }
649  double upCost = CoinMax((above-value)*upDynamicPseudoCost_,0.0);
650  return upCost;
651}
652// Return "down" estimate
653double 
654CbcSimpleIntegerDynamicPseudoCost::downEstimate() const
655{
656  const double * solution = model_->testSolution();
657  const double * lower = model_->getCbcColLower();
658  const double * upper = model_->getCbcColUpper();
659  double value = solution[columnNumber_];
660  value = CoinMax(value, lower[columnNumber_]);
661  value = CoinMin(value, upper[columnNumber_]);
662  if (upper[columnNumber_]==lower[columnNumber_]) {
663    // fixed
664    return 0.0;
665  }
666  double integerTolerance = 
667    model_->getDblParam(CbcModel::CbcIntegerTolerance);
668  double below = floor(value+integerTolerance);
669  double above = below+1.0;
670  if (above>upper[columnNumber_]) {
671    above=below;
672    below = above -1;
673  }
674  double downCost = CoinMax((value-below)*downDynamicPseudoCost_,0.0);
675  return downCost;
676}
677/* Pass in information on branch just done and create CbcObjectUpdateData instance.
678   If object does not need data then backward pointer will be NULL.
679   Assumes can get information from solver */
680CbcObjectUpdateData
681CbcSimpleIntegerDynamicPseudoCost::createUpdateInformation(const OsiSolverInterface * solver, 
682                                                           const CbcNode * node,
683                                                           const CbcBranchingObject * branchingObject)
684{
685  double originalValue=node->objectiveValue();
686  int originalUnsatisfied = node->numberUnsatisfied();
687  double objectiveValue = solver->getObjValue()*solver->getObjSense();
688  int unsatisfied=0;
689  int i;
690  //might be base model - doesn't matter
691  int numberIntegers = model_->numberIntegers();;
692  const double * solution = solver->getColSolution();
693  //const double * lower = solver->getColLower();
694  //const double * upper = solver->getColUpper();
695  double change = CoinMax(0.0,objectiveValue-originalValue);
696  int iStatus;
697  if (solver->isProvenOptimal())
698    iStatus=0; // optimal
699  else if (solver->isIterationLimitReached()
700           &&!solver->isDualObjectiveLimitReached())
701    iStatus=2; // unknown
702  else
703    iStatus=1; // infeasible
704
705  bool feasible = iStatus!=1;
706  if (feasible) {
707    double integerTolerance = 
708      model_->getDblParam(CbcModel::CbcIntegerTolerance);
709    const int * integerVariable = model_->integerVariable();
710    for (i=0;i<numberIntegers;i++) {
711      int j=integerVariable[i];
712      double value = solution[j];
713      double nearest = floor(value+0.5);
714      if (fabs(value-nearest)>integerTolerance) 
715        unsatisfied++;
716    }
717  }
718  int way = branchingObject->way();
719  way = - way; // because after branch so moved on
720  double value = branchingObject->value();
721  return CbcObjectUpdateData (this, way,
722                                  change, iStatus,
723                                  originalUnsatisfied-unsatisfied,value);
724}
725// Update object by CbcObjectUpdateData
726void 
727CbcSimpleIntegerDynamicPseudoCost::updateInformation(const CbcObjectUpdateData & data)
728{
729  bool feasible = data.status_!=1;
730  int way = data.way_;
731  double value = data.branchingValue_;
732  double change = data.change_;
733#define TYPE2 1
734#define TYPERATIO 0.9
735  if (way<0) {
736    // down
737    if (feasible) {
738      //printf("(down change %g value down %g ",change,value-floor(value));
739      incrementNumberTimesDown();
740      addToSumDownChange(1.0e-30+value-floor(value));
741      addToSumDownDecrease(data.intDecrease_);
742#if TYPE2==0
743      addToSumDownCost(change/(1.0e-30+(value-floor(value))));
744      setDownDynamicPseudoCost(sumDownCost()/(double) numberTimesDown());
745#elif TYPE2==1
746      addToSumDownCost(change);
747      setDownDynamicPseudoCost(sumDownCost()/sumDownChange());
748#elif TYPE2==2
749      addToSumDownCost(change*TYPERATIO+(1.0-TYPERATIO)*change/(1.0e-30+(value-floor(value))));
750      setDownDynamicPseudoCost(sumDownCost()*(TYPERATIO/sumDownChange()+(1.0-TYPERATIO)/(double) numberTimesDown()));
751#endif
752    } else {
753      //printf("(down infeasible value down %g ",change,value-floor(value));
754      incrementNumberTimesDownInfeasible();
755#if INFEAS==2
756      double distanceToCutoff=0.0;
757      double objectiveValue = model->getCurrentMinimizationObjValue();
758      distanceToCutoff =  model->getCutoff()  - originalValue;
759      if (distanceToCutoff<1.0e20) 
760        change = distanceToCutoff*2.0;
761      else 
762        change = downDynamicPseudoCost()*(value-floor(value))*10.0; 
763      change = CoinMax(1.0e-12*(1.0+fabs(originalValue)),change);
764      incrementNumberTimesDown();
765      addToSumDownChange(1.0e-30+value-floor(value));
766      addToSumDownDecrease(data.intDecrease_);
767#if TYPE2==0
768      addToSumDownCost(change/(1.0e-30+(value-floor(value))));
769      setDownDynamicPseudoCost(sumDownCost()/(double) numberTimesDown());
770#elif TYPE2==1
771      addToSumDownCost(change);
772      setDownDynamicPseudoCost(sumDownCost()/sumDownChange());
773#elif TYPE2==2
774      addToSumDownCost(change*TYPERATIO+(1.0-TYPERATIO)*change/(1.0e-30+(value-floor(value))));
775      setDownDynamicPseudoCost(sumDownCost()*(TYPERATIO/sumDownChange()+(1.0-TYPERATIO)/(double) numberTimesDown()));
776#endif
777#endif
778    }
779  } else {
780    // up
781    if (feasible) {
782      //printf("(up change %g value down %g ",change,ceil(value)-value);
783      incrementNumberTimesUp();
784      addToSumUpChange(1.0e-30+ceil(value)-value);
785      addToSumUpDecrease(data.intDecrease_);
786#if TYPE2==0
787      addToSumUpCost(change/(1.0e-30+(ceil(value)-value)));
788      setUpDynamicPseudoCost(sumUpCost()/(double) numberTimesUp());
789#elif TYPE2==1
790      addToSumUpCost(change);
791      setUpDynamicPseudoCost(sumUpCost()/sumUpChange());
792#elif TYPE2==2
793      addToSumUpCost(change*TYPERATIO+(1.0-TYPERATIO)*change/(1.0e-30+(ceil(value)-value)));
794      setUpDynamicPseudoCost(sumUpCost()*(TYPERATIO/sumUpChange()+(1.0-TYPERATIO)/(double) numberTimesUp()));
795#endif
796    } else {
797      //printf("(up infeasible value down %g ",change,ceil(value)-value);
798      incrementNumberTimesUpInfeasible();
799#if INFEAS==2
800      double distanceToCutoff=0.0;
801      double objectiveValue = model->getCurrentMinimizationObjValue();
802      distanceToCutoff =  model->getCutoff()  - originalValue;
803      if (distanceToCutoff<1.0e20) 
804        change = distanceToCutoff*2.0;
805      else 
806        change = upDynamicPseudoCost()*(ceil(value)-value)*10.0; 
807      change = CoinMax(1.0e-12*(1.0+fabs(originalValue)),change);
808      incrementNumberTimesUp();
809      addToSumUpChange(1.0e-30+ceil(value)-value);
810      addToSumUpDecrease(data.intDecrease_);
811#if TYPE2==0
812      addToSumUpCost(change/(1.0e-30+(ceil(value)-value)));
813      setUpDynamicPseudoCost(sumUpCost()/(double) numberTimesUp());
814#elif TYPE2==1
815      addToSumUpCost(change);
816      setUpDynamicPseudoCost(sumUpCost()/sumUpChange());
817#elif TYPE2==2
818      addToSumUpCost(change*TYPERATIO+(1.0-TYPERATIO)*change/(1.0e-30+(ceil(value)-value)));
819      setUpDynamicPseudoCost(sumUpCost()*(TYPERATIO/sumUpChange()+(1.0-TYPERATIO)/(double) numberTimesUp()));
820#endif
821#endif
822    }
823  }
824  //print(1,0.5);
825}
826// Pass in probing information
827void 
828CbcSimpleIntegerDynamicPseudoCost::setProbingInformation(int fixedDown, int fixedUp)
829{
830  numberTimesProbingTotal_++;
831  numberTimesDownLocalFixed_ = fixedDown;
832  numberTimesDownTotalFixed_ += fixedDown;
833  numberTimesUpLocalFixed_ = fixedUp;
834  numberTimesUpTotalFixed_ += fixedUp;
835}
836// Print
837void 
838CbcSimpleIntegerDynamicPseudoCost::print(int type,double value) const
839{
840  if (!type) {
841    double meanDown =0.0;
842    double devDown =0.0;
843    if (numberTimesDown_) {
844      meanDown = sumDownCost_/(double) numberTimesDown_;
845      devDown = meanDown*meanDown + sumDownCostSquared_ - 
846        2.0*meanDown*sumDownCost_;
847      if (devDown>=0.0)
848        devDown = sqrt(devDown);
849    }
850    double meanUp =0.0;
851    double devUp =0.0;
852    if (numberTimesUp_) {
853      meanUp = sumUpCost_/(double) numberTimesUp_;
854      devUp = meanUp*meanUp + sumUpCostSquared_ - 
855        2.0*meanUp*sumUpCost_;
856      if (devUp>=0.0)
857        devUp = sqrt(devUp);
858    }
859#if 0
860    printf("%d down %d times (%d inf) mean %g (dev %g) up %d times (%d inf) mean %g (dev %g)\n",
861           columnNumber_,
862           numberTimesDown_,numberTimesDownInfeasible_,meanDown,devDown,
863           numberTimesUp_,numberTimesUpInfeasible_,meanUp,devUp);
864#else
865    printf("%d down %d times (%d inf) mean %g  up %d times (%d inf) mean %g\n",
866           columnNumber_,
867           numberTimesDown_,numberTimesDownInfeasible_,meanDown,
868           numberTimesUp_,numberTimesUpInfeasible_,meanUp);
869#endif
870  } else {
871    const double * upper = model_->getCbcColUpper();
872    double integerTolerance = 
873      model_->getDblParam(CbcModel::CbcIntegerTolerance);
874    double below = floor(value+integerTolerance);
875    double above = below+1.0;
876    if (above>upper[columnNumber_]) {
877      above=below;
878      below = above -1;
879    }
880    double objectiveValue = model_->getCurrentMinimizationObjValue();
881    double distanceToCutoff =  model_->getCutoff() - objectiveValue;
882    if (distanceToCutoff<1.0e20) 
883      distanceToCutoff *= 10.0;
884    else 
885      distanceToCutoff = 1.0e2 + fabs(objectiveValue);
886    double sum;
887    int number;
888    double downCost = CoinMax(value-below,0.0);
889    double downCost0 = downCost*downDynamicPseudoCost_;
890    sum = sumDownCost();
891    number = numberTimesDown();
892    sum += numberTimesDownInfeasible()*(distanceToCutoff/(downCost+1.0e-12));
893    if (number>0)
894      downCost *= sum / (double) number;
895    else
896      downCost  *=  downDynamicPseudoCost_;
897    double upCost = CoinMax((above-value),0.0);
898    double upCost0 = upCost*upDynamicPseudoCost_;
899    sum = sumUpCost();
900    number = numberTimesUp();
901    sum += numberTimesUpInfeasible()*(distanceToCutoff/(upCost+1.0e-12));
902    if (number>0)
903      upCost *= sum / (double) number;
904    else
905      upCost  *=  upDynamicPseudoCost_;
906    printf("%d down %d times %g (est %g)  up %d times %g (est %g)\n",
907           columnNumber_,
908           numberTimesDown_,downCost,downCost0,
909           numberTimesUp_,upCost,upCost0);
910  }
911}
912
913// Default Constructor
914CbcDynamicPseudoCostBranchingObject::CbcDynamicPseudoCostBranchingObject()
915  :CbcIntegerBranchingObject()
916{
917  changeInGuessed_=1.0e-5;
918  object_=NULL;
919}
920
921// Useful constructor
922CbcDynamicPseudoCostBranchingObject::CbcDynamicPseudoCostBranchingObject (CbcModel * model, 
923                                                                          int variable, 
924                                                                          int way , double value,
925                                       CbcSimpleIntegerDynamicPseudoCost * object) 
926  :CbcIntegerBranchingObject(model,variable,way,value)
927{
928  changeInGuessed_=1.0e-5;
929  object_=object;
930}
931// Useful constructor for fixing
932CbcDynamicPseudoCostBranchingObject::CbcDynamicPseudoCostBranchingObject (CbcModel * model, 
933                                                      int variable, int way,
934                                                      double lowerValue, 
935                                                      double upperValue)
936  :CbcIntegerBranchingObject(model,variable,way,lowerValue)
937{
938  changeInGuessed_=1.0e100;
939  object_=NULL;
940}
941 
942
943// Copy constructor
944CbcDynamicPseudoCostBranchingObject::CbcDynamicPseudoCostBranchingObject ( 
945                                 const CbcDynamicPseudoCostBranchingObject & rhs)
946  :CbcIntegerBranchingObject(rhs)
947{
948  changeInGuessed_ = rhs.changeInGuessed_;
949  object_=rhs.object_;
950}
951
952// Assignment operator
953CbcDynamicPseudoCostBranchingObject & 
954CbcDynamicPseudoCostBranchingObject::operator=( const CbcDynamicPseudoCostBranchingObject& rhs)
955{
956  if (this != &rhs) {
957    CbcIntegerBranchingObject::operator=(rhs);
958    changeInGuessed_ = rhs.changeInGuessed_;
959    object_=rhs.object_;
960  }
961  return *this;
962}
963CbcBranchingObject * 
964CbcDynamicPseudoCostBranchingObject::clone() const
965{ 
966  return (new CbcDynamicPseudoCostBranchingObject(*this));
967}
968
969
970// Destructor
971CbcDynamicPseudoCostBranchingObject::~CbcDynamicPseudoCostBranchingObject ()
972{
973}
974
975/*
976  Perform a branch by adjusting the bounds of the specified variable. Note
977  that each arm of the branch advances the object to the next arm by
978  advancing the value of way_.
979
980  Providing new values for the variable's lower and upper bounds for each
981  branching direction gives a little bit of additional flexibility and will
982  be easily extensible to multi-way branching.
983  Returns change in guessed objective on next branch
984*/
985double
986CbcDynamicPseudoCostBranchingObject::branch()
987{
988  CbcIntegerBranchingObject::branch();
989  return changeInGuessed_;
990}
991/* Some branchingObjects may claim to be able to skip
992   strong branching.  If so they have to fill in CbcStrongInfo.
993   The object mention in incoming CbcStrongInfo must match.
994   Returns nonzero if skip is wanted */
995int 
996CbcDynamicPseudoCostBranchingObject::fillStrongInfo( CbcStrongInfo & info)
997{
998  assert (object_);
999  assert (info.possibleBranch==this);
1000  if (object_->numberTimesUp()<object_->numberBeforeTrust()||
1001      object_->numberTimesDown()<object_->numberBeforeTrust()) {
1002    return 0;
1003  } else {
1004    info.upMovement = object_->upDynamicPseudoCost()*(ceil(value_)-value_);
1005    info.downMovement = object_->downDynamicPseudoCost()*(value_-floor(value_));
1006    info.numIntInfeasUp  -= (int) (object_->sumUpDecrease()/
1007                                   ((double) object_->numberTimesUp()));
1008    info.numObjInfeasUp = 0;
1009    info.finishedUp = false;
1010    info.numItersUp = 0;
1011    info.numIntInfeasDown  -= (int) (object_->sumDownDecrease()/
1012                                   ((double) object_->numberTimesDown()));
1013    info.numObjInfeasDown = 0;
1014    info.finishedDown = false;
1015    info.numItersDown = 0;
1016    info.fix =0;
1017    return 1;
1018  }
1019}
1020 
1021// Default Constructor
1022CbcBranchDynamicDecision::CbcBranchDynamicDecision()
1023  :CbcBranchDecision()
1024{
1025  bestCriterion_ = 0.0;
1026  bestChangeUp_ = 0.0;
1027  bestNumberUp_ = 0;
1028  bestChangeDown_ = 0.0;
1029  bestNumberDown_ = 0;
1030  bestObject_ = NULL;
1031}
1032
1033// Copy constructor
1034CbcBranchDynamicDecision::CbcBranchDynamicDecision (
1035                                    const CbcBranchDynamicDecision & rhs)
1036  :CbcBranchDecision()
1037{
1038  bestCriterion_ = rhs.bestCriterion_;
1039  bestChangeUp_ = rhs.bestChangeUp_;
1040  bestNumberUp_ = rhs.bestNumberUp_;
1041  bestChangeDown_ = rhs.bestChangeDown_;
1042  bestNumberDown_ = rhs.bestNumberDown_;
1043  bestObject_ = rhs.bestObject_;
1044}
1045
1046CbcBranchDynamicDecision::~CbcBranchDynamicDecision()
1047{
1048}
1049
1050// Clone
1051CbcBranchDecision * 
1052CbcBranchDynamicDecision::clone() const
1053{
1054  return new CbcBranchDynamicDecision(*this);
1055}
1056
1057// Initialize i.e. before start of choosing at a node
1058void 
1059CbcBranchDynamicDecision::initialize(CbcModel * model)
1060{
1061  bestCriterion_ = 0.0;
1062  bestChangeUp_ = 0.0;
1063  bestNumberUp_ = 0;
1064  bestChangeDown_ = 0.0;
1065  bestNumberDown_ = 0;
1066  bestObject_ = NULL;
1067}
1068
1069/* Saves a clone of current branching object.  Can be used to update
1070      information on object causing branch - after branch */
1071void 
1072CbcBranchDynamicDecision::saveBranchingObject(OsiBranchingObject * object) 
1073{
1074  OsiBranchingObject * obj = object->clone();
1075  CbcDynamicPseudoCostBranchingObject * branchingObject =
1076    dynamic_cast<CbcDynamicPseudoCostBranchingObject *>(obj);
1077  assert (branchingObject);
1078  object_=branchingObject;
1079}
1080/* Pass in information on branch just done.
1081   assumes object can get information from solver */
1082void 
1083CbcBranchDynamicDecision::updateInformation(OsiSolverInterface * solver,
1084                                            const CbcNode * node)
1085{
1086  assert (object_);
1087  const CbcModel * model = object_->model();
1088  double originalValue=node->objectiveValue();
1089  int originalUnsatisfied = node->numberUnsatisfied();
1090  double objectiveValue = solver->getObjValue()*model->getObjSense();
1091  int unsatisfied=0;
1092  int i;
1093  int numberIntegers = model->numberIntegers();;
1094  const double * solution = solver->getColSolution();
1095  //const double * lower = solver->getColLower();
1096  //const double * upper = solver->getColUpper();
1097  CbcDynamicPseudoCostBranchingObject * branchingObject =
1098    dynamic_cast<CbcDynamicPseudoCostBranchingObject *>(object_);
1099  if (!branchingObject) {
1100    delete object_;
1101    object_=NULL;
1102    return;
1103  }
1104  CbcSimpleIntegerDynamicPseudoCost *  object = branchingObject->object();
1105  double change = CoinMax(0.0,objectiveValue-originalValue);
1106  // probably should also ignore if stopped
1107  int iStatus;
1108  if (solver->isProvenOptimal())
1109    iStatus=0; // optimal
1110  else if (solver->isIterationLimitReached()
1111           &&!solver->isDualObjectiveLimitReached())
1112    iStatus=2; // unknown
1113  else
1114    iStatus=1; // infeasible
1115
1116  bool feasible = iStatus!=1;
1117  if (feasible) {
1118    double integerTolerance = 
1119      model->getDblParam(CbcModel::CbcIntegerTolerance);
1120    const int * integerVariable = model->integerVariable();
1121    for (i=0;i<numberIntegers;i++) {
1122      int j=integerVariable[i];
1123      double value = solution[j];
1124      double nearest = floor(value+0.5);
1125      if (fabs(value-nearest)>integerTolerance) 
1126        unsatisfied++;
1127    }
1128  }
1129  int way = object_->way();
1130  double value = object_->value();
1131  //#define TYPE2 1
1132  //#define TYPERATIO 0.9
1133  if (way<0) {
1134    // down
1135    if (feasible) {
1136      //printf("(down change %g value down %g ",change,value-floor(value));
1137      object->incrementNumberTimesDown();
1138      object->addToSumDownChange(1.0e-30+value-floor(value));
1139      object->addToSumDownDecrease(originalUnsatisfied-unsatisfied);
1140#if TYPE2==0
1141      object->addToSumDownCost(change/(1.0e-30+(value-floor(value))));
1142      object->setDownDynamicPseudoCost(object->sumDownCost()/(double) object->numberTimesDown());
1143#elif TYPE2==1
1144      object->addToSumDownCost(change);
1145      object->setDownDynamicPseudoCost(object->sumDownCost()/object->sumDownChange());
1146#elif TYPE2==2
1147      object->addToSumDownCost(change*TYPERATIO+(1.0-TYPERATIO)*change/(1.0e-30+(value-floor(value))));
1148      object->setDownDynamicPseudoCost(object->sumDownCost()*(TYPERATIO/object->sumDownChange()+(1.0-TYPERATIO)/(double) object->numberTimesDown()));
1149#endif
1150    } else {
1151      //printf("(down infeasible value down %g ",change,value-floor(value));
1152      object->incrementNumberTimesDownInfeasible();
1153#if INFEAS==2
1154      double distanceToCutoff=0.0;
1155      double objectiveValue = model->getCurrentMinimizationObjValue();
1156      distanceToCutoff =  model->getCutoff()  - originalValue;
1157      if (distanceToCutoff<1.0e20) 
1158        change = distanceToCutoff*2.0;
1159      else 
1160        change = object->downDynamicPseudoCost()*(value-floor(value))*10.0; 
1161      change = CoinMax(1.0e-12*(1.0+fabs(originalValue)),change);
1162      object->incrementNumberTimesDown();
1163      object->addToSumDownChange(1.0e-30+value-floor(value));
1164      object->addToSumDownDecrease(originalUnsatisfied-unsatisfied);
1165#if TYPE2==0
1166      object->addToSumDownCost(change/(1.0e-30+(value-floor(value))));
1167      object->setDownDynamicPseudoCost(object->sumDownCost()/(double) object->numberTimesDown());
1168#elif TYPE2==1
1169      object->addToSumDownCost(change);
1170      object->setDownDynamicPseudoCost(object->sumDownCost()/object->sumDownChange());
1171#elif TYPE2==2
1172      object->addToSumDownCost(change*TYPERATIO+(1.0-TYPERATIO)*change/(1.0e-30+(value-floor(value))));
1173      object->setDownDynamicPseudoCost(object->sumDownCost()*(TYPERATIO/object->sumDownChange()+(1.0-TYPERATIO)/(double) object->numberTimesDown()));
1174#endif
1175#endif
1176    }
1177  } else {
1178    // up
1179    if (feasible) {
1180      //printf("(up change %g value down %g ",change,ceil(value)-value);
1181      object->incrementNumberTimesUp();
1182      object->addToSumUpChange(1.0e-30+ceil(value)-value);
1183      object->addToSumUpDecrease(unsatisfied-originalUnsatisfied);
1184#if TYPE2==0
1185      object->addToSumUpCost(change/(1.0e-30+(ceil(value)-value)));
1186      object->setUpDynamicPseudoCost(object->sumUpCost()/(double) object->numberTimesUp());
1187#elif TYPE2==1
1188      object->addToSumUpCost(change);
1189      object->setUpDynamicPseudoCost(object->sumUpCost()/object->sumUpChange());
1190#elif TYPE2==2
1191      object->addToSumUpCost(change*TYPERATIO+(1.0-TYPERATIO)*change/(1.0e-30+(ceil(value)-value)));
1192      object->setUpDynamicPseudoCost(object->sumUpCost()*(TYPERATIO/object->sumUpChange()+(1.0-TYPERATIO)/(double) object->numberTimesUp()));
1193#endif
1194    } else {
1195      //printf("(up infeasible value down %g ",change,ceil(value)-value);
1196      object->incrementNumberTimesUpInfeasible();
1197#if INFEAS==2
1198      double distanceToCutoff=0.0;
1199      double objectiveValue = model->getCurrentMinimizationObjValue();
1200      distanceToCutoff =  model->getCutoff()  - originalValue;
1201      if (distanceToCutoff<1.0e20) 
1202        change = distanceToCutoff*2.0;
1203      else 
1204        change = object->upDynamicPseudoCost()*(ceil(value)-value)*10.0; 
1205      change = CoinMax(1.0e-12*(1.0+fabs(originalValue)),change);
1206      object->incrementNumberTimesUp();
1207      object->addToSumUpChange(1.0e-30+ceil(value)-value);
1208      object->addToSumUpDecrease(unsatisfied-originalUnsatisfied);
1209#if TYPE2==0
1210      object->addToSumUpCost(change/(1.0e-30+(ceil(value)-value)));
1211      object->setUpDynamicPseudoCost(object->sumUpCost()/(double) object->numberTimesUp());
1212#elif TYPE2==1
1213      object->addToSumUpCost(change);
1214      object->setUpDynamicPseudoCost(object->sumUpCost()/object->sumUpChange());
1215#elif TYPE2==2
1216      object->addToSumUpCost(change*TYPERATIO+(1.0-TYPERATIO)*change/(1.0e-30+(ceil(value)-value)));
1217      object->setUpDynamicPseudoCost(object->sumUpCost()*(TYPERATIO/object->sumUpChange()+(1.0-TYPERATIO)/(double) object->numberTimesUp()));
1218#endif
1219#endif
1220    }
1221  }
1222  //object->print(1,0.5);
1223  delete object_;
1224  object_=NULL;
1225}
1226
1227/*
1228  Simple dynamic decision algorithm. Compare based on infeasibility (numInfUp,
1229  numInfDown) until a solution is found by search, then switch to change in
1230  objective (changeUp, changeDown). Note that bestSoFar is remembered in
1231  bestObject_, so the parameter bestSoFar is unused.
1232*/
1233
1234int
1235CbcBranchDynamicDecision::betterBranch(CbcBranchingObject * thisOne,
1236                            CbcBranchingObject * bestSoFar,
1237                            double changeUp, int numInfUp,
1238                            double changeDown, int numInfDown)
1239{
1240  CbcModel * model = thisOne->model();
1241  int stateOfSearch = model->stateOfSearch();
1242  int betterWay=0;
1243  double value=0.0;
1244  if (!bestObject_) {
1245    bestCriterion_=-1.0;
1246    bestNumberUp_=INT_MAX;
1247    bestNumberDown_=INT_MAX;
1248  }
1249  if (stateOfSearch<=1&&thisOne->model()->currentNode()->depth()>10) {
1250#define TRY_STUFF 1
1251#ifdef TRY_STUFF
1252    // before solution - choose smallest number
1253    // could add in depth as well
1254    int bestNumber = CoinMin(bestNumberUp_,bestNumberDown_);
1255    if (numInfUp<numInfDown) {
1256      if (numInfUp<bestNumber) {
1257        betterWay = 1;
1258      } else if (numInfUp==bestNumber) {
1259        if (changeUp<bestChangeUp_)
1260          betterWay=1;
1261      }
1262    } else if (numInfUp>numInfDown) {
1263      if (numInfDown<bestNumber) {
1264        betterWay = -1;
1265      } else if (numInfDown==bestNumber) {
1266        if (changeDown<bestChangeDown_)
1267          betterWay=-1;
1268      }
1269    } else {
1270      // up and down have same number
1271      bool better=false;
1272      if (numInfUp<bestNumber) {
1273        better=true;
1274      } else if (numInfUp==bestNumber) {
1275        if (CoinMin(changeUp,changeDown)<CoinMin(bestChangeUp_,bestChangeDown_)-1.0e-5)
1276          better=true;;
1277      }
1278      if (better) {
1279        // see which way
1280        if (changeUp<=changeDown)
1281          betterWay=1;
1282        else
1283          betterWay=-1;
1284      }
1285    }
1286    if (betterWay) {
1287      value = CoinMin(numInfUp,numInfDown);
1288    }
1289#else
1290    // got a solution
1291    double minValue = CoinMin(changeDown,changeUp);
1292    double maxValue = CoinMax(changeDown,changeUp);
1293    value = WEIGHT_BEFORE*minValue + (1.0-WEIGHT_BEFORE)*maxValue;
1294    if (value>bestCriterion_+1.0e-8) {
1295      if (changeUp<=changeDown) {
1296        betterWay=1;
1297      } else {
1298        betterWay=-1;
1299      }
1300    }
1301#endif
1302  } else {
1303#if TRY_STUFF > 1
1304    // Get current number of infeasibilities, cutoff and current objective
1305    CbcNode * node = model->currentNode();
1306    int numberUnsatisfied = node->numberUnsatisfied();
1307    double cutoff = model->getCutoff();
1308    double objectiveValue = node->objectiveValue();
1309#endif
1310    // got a solution
1311    double minValue = CoinMin(changeDown,changeUp);
1312    double maxValue = CoinMax(changeDown,changeUp);
1313    // Reduce
1314#ifdef TRY_STUFF
1315    //maxValue = CoinMin(maxValue,minValue*4.0);
1316#else
1317    maxValue = CoinMin(maxValue,minValue*2.0);
1318#endif
1319    value = WEIGHT_AFTER*minValue + (1.0-WEIGHT_AFTER)*maxValue;
1320    double useValue = value;
1321    double useBest = bestCriterion_;
1322#if TRY_STUFF>1
1323    int thisNumber = CoinMin(numInfUp,numInfDown);
1324    int bestNumber = CoinMin(bestNumberUp_,bestNumberDown_);
1325    double distance = cutoff-objectiveValue;
1326    assert (distance>=0.0);
1327    if (useValue+0.1*distance>useBest&&useValue*1.1>useBest&&
1328        useBest+0.1*distance>useValue&&useBest*1.1>useValue) {
1329      // not much in it - look at unsatisfied
1330      if (thisNumber<numberUnsatisfied||bestNumber<numberUnsatisfied) {
1331        double perInteger = distance/ ((double) numberUnsatisfied);
1332        useValue += thisNumber*perInteger;
1333        useBest += bestNumber*perInteger;
1334      }
1335    }
1336#endif
1337    if (useValue>useBest+1.0e-8) {
1338      if (changeUp<=changeDown) {
1339        betterWay=1;
1340      } else {
1341        betterWay=-1;
1342      }
1343    }
1344  }
1345  if (betterWay) {
1346    // maybe change better way
1347    CbcDynamicPseudoCostBranchingObject * branchingObject =
1348      dynamic_cast<CbcDynamicPseudoCostBranchingObject *>(thisOne);
1349    if (branchingObject) {
1350      CbcSimpleIntegerDynamicPseudoCost *  object = branchingObject->object();
1351      double separator = object->upDownSeparator();
1352      if (separator>0.0) {
1353        const double * solution = thisOne->model()->testSolution();
1354        double valueVariable = solution[object->columnNumber()];
1355        betterWay = (valueVariable-floor(valueVariable)>=separator) ? 1 : -1;
1356      }
1357    }
1358    bestCriterion_ = value;
1359    bestChangeUp_ = changeUp;
1360    bestNumberUp_ = numInfUp;
1361    bestChangeDown_ = changeDown;
1362    bestNumberDown_ = numInfDown;
1363    bestObject_=thisOne;
1364    // See if user is overriding way
1365    if (thisOne->object()&&thisOne->object()->preferredWay())
1366      betterWay = thisOne->object()->preferredWay();
1367  }
1368  return betterWay;
1369}
1370/* Sets or gets best criterion so far */
1371void 
1372CbcBranchDynamicDecision::setBestCriterion(double value)
1373{ 
1374  bestCriterion_ = value;
1375}
1376double 
1377CbcBranchDynamicDecision::getBestCriterion() const
1378{ 
1379  return bestCriterion_;
1380}
Note: See TracBrowser for help on using the repository browser.