source: stable/2.0/Cbc/src/CbcBranchDynamic.cpp @ 905

Last change on this file since 905 was 905, checked in by ladanyi, 13 years ago

include cstdlib before cmath to get things to compile on AIX with xlC (same as changeset 904 in trunk)

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