source: trunk/Cbc/src/CbcBranchDynamic.cpp @ 702

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

changes for gcc 4.3

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 46.7 KB
Line 
1// Copyright (C) 2002, International Business Machines
2// Corporation and others.  All Rights Reserved.
3#if defined(_MSC_VER)
4// Turn off compiler warning about long names
5#  pragma warning(disable:4786)
6#endif
7#include <cassert>
8#include <cmath>
9#include <cfloat>
10//#define CBC_DEBUG
11
12#include "OsiSolverInterface.hpp"
13#include "OsiSolverBranch.hpp"
14#include "CbcModel.hpp"
15#include "CbcMessage.hpp"
16#include "CbcBranchDynamic.hpp"
17#include "CoinSort.hpp"
18#include "CoinError.hpp"
19
20/** Default Constructor
21
22  Equivalent to an unspecified binary variable.
23*/
24CbcSimpleIntegerDynamicPseudoCost::CbcSimpleIntegerDynamicPseudoCost ()
25  : CbcSimpleInteger(),
26    downDynamicPseudoCost_(1.0e-5),
27    upDynamicPseudoCost_(1.0e-5),
28    upDownSeparator_(-1.0),
29    sumDownCost_(0.0),
30    sumUpCost_(0.0),
31    sumDownChange_(0.0),
32    sumUpChange_(0.0),
33    sumDownCostSquared_(0.0),
34    sumUpCostSquared_(0.0),
35    sumDownDecrease_(0.0),
36    sumUpDecrease_(0.0),
37    lastDownCost_(0.0),
38    lastUpCost_(0.0),
39    lastDownDecrease_(0),
40    lastUpDecrease_(0),
41    numberTimesDown_(0),
42    numberTimesUp_(0),
43    numberTimesDownInfeasible_(0),
44    numberTimesUpInfeasible_(0),
45    numberBeforeTrust_(0),
46    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.7
437#define WEIGHT_BEFORE 0.1
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      if (0) {
449        double sum;
450        int number;
451        double downCost2 = CoinMax(value-below,0.0);
452        sum = sumDownCost_;
453        number = numberTimesDown_;
454        if (number>0)
455          downCost2 *= sum / (double) number;
456        else
457          downCost2  *=  downDynamicPseudoCost_;
458        double upCost2 = CoinMax((above-value),0.0);
459        sum = sumUpCost_;
460        number = numberTimesUp_;
461        if (number>0)
462          upCost2 *= sum / (double) number;
463        else
464          upCost2  *=  upDynamicPseudoCost_;
465        double minValue2 = CoinMin(downCost2,upCost2);
466        double maxValue2 = CoinMax(downCost2,upCost2);
467        printf("%d value %g downC %g upC %g minV %g maxV %g downC2 %g upC2 %g minV2 %g maxV2 %g\n",
468               columnNumber_,value,downCost,upCost,minValue,maxValue,
469               downCost2,upCost2,minValue2,maxValue2);
470      }
471    } else {
472      // some solution
473      returnValue = WEIGHT_AFTER*minValue + (1.0-WEIGHT_AFTER)*maxValue;
474    }
475    if (numberTimesUp_<numberBeforeTrust_||
476        numberTimesDown_<numberBeforeTrust_) {
477      //if (returnValue<1.0e10)
478      //returnValue += 1.0e12;
479      //else
480      returnValue *= 1.0e3;
481      if (!numberTimesUp_&&!numberTimesDown_)
482        returnValue=1.0e50;
483    }
484    //if (fabs(value-0.5)<1.0e-5) {
485    //returnValue = 3.0*returnValue + 0.2;
486    //} else if (value>0.9) {
487    //returnValue = 2.0*returnValue + 0.1;
488    //}
489    if (method_==1) {
490      // probing
491      // average
492      double up=1.0e-15;
493      double down=1.0e-15;
494      if (numberTimesProbingTotal_) {
495        up += numberTimesUpTotalFixed_/((double) numberTimesProbingTotal_);
496        down += numberTimesDownTotalFixed_/((double) numberTimesProbingTotal_);
497      }
498      returnValue = 1 + 10.0*CoinMin(numberTimesDownLocalFixed_,numberTimesUpLocalFixed_) +
499        CoinMin(down,up);
500      returnValue *= 1.0e-3;
501    }
502    return CoinMax(returnValue,1.0e-15);
503  }
504}
505
506double
507CbcSimpleIntegerDynamicPseudoCost::infeasibility(const OsiSolverInterface * solver, const OsiBranchingInformation * info,
508                         int & preferredWay) const
509{
510  double value = info->solution_[columnNumber_];
511  value = CoinMax(value, info->lower_[columnNumber_]);
512  value = CoinMin(value, info->upper_[columnNumber_]);
513  if (info->upper_[columnNumber_]==info->lower_[columnNumber_]) {
514    // fixed
515    preferredWay=1;
516    return 0.0;
517  }
518  assert (breakEven_>0.0&&breakEven_<1.0);
519  double nearest = floor(value+0.5);
520  double integerTolerance = info->integerTolerance_; 
521  double below = floor(value+integerTolerance);
522  double above = below+1.0;
523  if (above>info->upper_[columnNumber_]) {
524    above=below;
525    below = above -1;
526  }
527#if INFEAS==1
528  double objectiveValue = info->objectiveValue_;
529  double distanceToCutoff =  info->cutoff_  - objectiveValue;
530  if (distanceToCutoff<1.0e20) 
531    distanceToCutoff *= 10.0;
532  else 
533    distanceToCutoff = 1.0e2 + fabs(objectiveValue);
534#endif
535  double sum;
536  int number;
537  double downCost = CoinMax(value-below,0.0);
538  sum = sumDownCost_;
539  number = numberTimesDown_;
540#if INFEAS==1
541  sum += numberTimesDownInfeasible_*(distanceToCutoff/(downCost+1.0e-12));
542  number += numberTimesDownInfeasible_;
543#endif
544  if (number>0)
545    downCost *= sum / (double) number;
546  else
547    downCost  *=  downDynamicPseudoCost_;
548  double upCost = CoinMax((above-value),0.0);
549  sum = sumUpCost_;
550  number = numberTimesUp_;
551#if INFEAS==1
552  sum += numberTimesUpInfeasible_*(distanceToCutoff/(upCost+1.0e-12));
553  number += numberTimesUpInfeasible_;
554#endif
555  if (number>0)
556    upCost *= sum / (double) number;
557  else
558    upCost  *=  upDynamicPseudoCost_;
559  if (downCost>=upCost)
560    preferredWay=1;
561  else
562    preferredWay=-1;
563  // See if up down choice set
564  if (upDownSeparator_>0.0) {
565    preferredWay = (value-below>=upDownSeparator_) ? 1 : -1;
566  }
567  if (preferredWay_)
568    preferredWay=preferredWay_;
569  // weight at 1.0 is max min
570  if (fabs(value-nearest)<=integerTolerance) {
571    return 0.0;
572  } else {
573    double returnValue=0.0;
574    double minValue = CoinMin(downCost,upCost);
575    double maxValue = CoinMax(downCost,upCost);
576    if (!info->numberBranchingSolutions_||info->depth_<=10/* was ||maxValue>0.2*distanceToCutoff*/) {
577      // no solution
578      returnValue = WEIGHT_BEFORE*minValue + (1.0-WEIGHT_BEFORE)*maxValue;
579    } else {
580      // some solution
581      returnValue = WEIGHT_AFTER*minValue + (1.0-WEIGHT_AFTER)*maxValue;
582    }
583    if (numberTimesUp_<numberBeforeTrust_||
584        numberTimesDown_<numberBeforeTrust_) {
585      //if (returnValue<1.0e10)
586      //returnValue += 1.0e12;
587      //else
588      returnValue *= 1.0e3;
589      if (!numberTimesUp_&&!numberTimesDown_)
590        returnValue=1.0e50;
591    }
592    //if (fabs(value-0.5)<1.0e-5) {
593    //returnValue = 3.0*returnValue + 0.2;
594    //} else if (value>0.9) {
595    //returnValue = 2.0*returnValue + 0.1;
596    //}
597    if (method_==1) {
598      // probing
599      // average
600      double up=1.0e-15;
601      double down=1.0e-15;
602      if (numberTimesProbingTotal_) {
603        up += numberTimesUpTotalFixed_/((double) numberTimesProbingTotal_);
604        down += numberTimesDownTotalFixed_/((double) numberTimesProbingTotal_);
605      }
606      returnValue = 1 + 10.0*CoinMin(numberTimesDownLocalFixed_,numberTimesUpLocalFixed_) +
607        CoinMin(down,up);
608      returnValue *= 1.0e-3;
609    }
610    return CoinMax(returnValue,1.0e-15);
611  }
612}
613// Creates a branching object
614CbcBranchingObject * 
615CbcSimpleIntegerDynamicPseudoCost::createBranch(OsiSolverInterface * solver, const OsiBranchingInformation * info, int way) 
616{
617  double value = info->solution_[columnNumber_];
618  value = CoinMax(value, info->lower_[columnNumber_]);
619  value = CoinMin(value, info->upper_[columnNumber_]);
620  assert (info->upper_[columnNumber_]>info->lower_[columnNumber_]);
621  if (!info->hotstartSolution_) {
622#ifndef NDEBUG
623    double nearest = floor(value+0.5);
624    assert (fabs(value-nearest)>info->integerTolerance_);
625#endif
626  } else {
627    double targetValue = info->hotstartSolution_[columnNumber_];
628    if (way>0)
629      value = targetValue-0.1;
630    else
631      value = targetValue+0.1;
632  }
633  CbcDynamicPseudoCostBranchingObject * newObject = 
634    new CbcDynamicPseudoCostBranchingObject(model_,columnNumber_,way,
635                                            value,this);
636  double up =  upDynamicPseudoCost_*(ceil(value)-value);
637  double down =  downDynamicPseudoCost_*(value-floor(value));
638  double changeInGuessed=up-down;
639  if (way>0)
640    changeInGuessed = - changeInGuessed;
641  changeInGuessed=CoinMax(0.0,changeInGuessed);
642  //if (way>0)
643  //changeInGuessed += 1.0e8; // bias to stay up
644  newObject->setChangeInGuessed(changeInGuessed);
645  newObject->setOriginalObject(this);
646  return newObject;
647}
648
649// Return "up" estimate
650double 
651CbcSimpleIntegerDynamicPseudoCost::upEstimate() const
652{
653  const double * solution = model_->testSolution();
654  const double * lower = model_->getCbcColLower();
655  const double * upper = model_->getCbcColUpper();
656  double value = solution[columnNumber_];
657  value = CoinMax(value, lower[columnNumber_]);
658  value = CoinMin(value, upper[columnNumber_]);
659  if (upper[columnNumber_]==lower[columnNumber_]) {
660    // fixed
661    return 0.0;
662  }
663  double integerTolerance = 
664    model_->getDblParam(CbcModel::CbcIntegerTolerance);
665  double below = floor(value+integerTolerance);
666  double above = below+1.0;
667  if (above>upper[columnNumber_]) {
668    above=below;
669    below = above -1;
670  }
671  double upCost = CoinMax((above-value)*upDynamicPseudoCost_,0.0);
672  return upCost;
673}
674// Return "down" estimate
675double 
676CbcSimpleIntegerDynamicPseudoCost::downEstimate() const
677{
678  const double * solution = model_->testSolution();
679  const double * lower = model_->getCbcColLower();
680  const double * upper = model_->getCbcColUpper();
681  double value = solution[columnNumber_];
682  value = CoinMax(value, lower[columnNumber_]);
683  value = CoinMin(value, upper[columnNumber_]);
684  if (upper[columnNumber_]==lower[columnNumber_]) {
685    // fixed
686    return 0.0;
687  }
688  double integerTolerance = 
689    model_->getDblParam(CbcModel::CbcIntegerTolerance);
690  double below = floor(value+integerTolerance);
691  double above = below+1.0;
692  if (above>upper[columnNumber_]) {
693    above=below;
694    below = above -1;
695  }
696  double downCost = CoinMax((value-below)*downDynamicPseudoCost_,0.0);
697  return downCost;
698}
699/* Pass in information on branch just done and create CbcObjectUpdateData instance.
700   If object does not need data then backward pointer will be NULL.
701   Assumes can get information from solver */
702CbcObjectUpdateData
703CbcSimpleIntegerDynamicPseudoCost::createUpdateInformation(const OsiSolverInterface * solver, 
704                                                           const CbcNode * node,
705                                                           const CbcBranchingObject * branchingObject)
706{
707  double originalValue=node->objectiveValue();
708  int originalUnsatisfied = node->numberUnsatisfied();
709  double objectiveValue = solver->getObjValue()*solver->getObjSense();
710  int unsatisfied=0;
711  int i;
712  //might be base model - doesn't matter
713  int numberIntegers = model_->numberIntegers();;
714  const double * solution = solver->getColSolution();
715  //const double * lower = solver->getColLower();
716  //const double * upper = solver->getColUpper();
717  double change = CoinMax(0.0,objectiveValue-originalValue);
718  int iStatus;
719  if (solver->isProvenOptimal())
720    iStatus=0; // optimal
721  else if (solver->isIterationLimitReached()
722           &&!solver->isDualObjectiveLimitReached())
723    iStatus=2; // unknown
724  else
725    iStatus=1; // infeasible
726
727  bool feasible = iStatus!=1;
728  if (feasible) {
729    double integerTolerance = 
730      model_->getDblParam(CbcModel::CbcIntegerTolerance);
731    const int * integerVariable = model_->integerVariable();
732    for (i=0;i<numberIntegers;i++) {
733      int j=integerVariable[i];
734      double value = solution[j];
735      double nearest = floor(value+0.5);
736      if (fabs(value-nearest)>integerTolerance) 
737        unsatisfied++;
738    }
739  }
740  int way = branchingObject->way();
741  way = - way; // because after branch so moved on
742  double value = branchingObject->value();
743  return CbcObjectUpdateData (this, way,
744                                  change, iStatus,
745                                  originalUnsatisfied-unsatisfied,value);
746}
747// Update object by CbcObjectUpdateData
748void 
749CbcSimpleIntegerDynamicPseudoCost::updateInformation(const CbcObjectUpdateData & data)
750{
751  bool feasible = data.status_!=1;
752  int way = data.way_;
753  double value = data.branchingValue_;
754  double change = data.change_;
755#define TYPE2 1
756#define TYPERATIO 0.9
757  if (way<0) {
758    // down
759    if (feasible) {
760      //printf("(down change %g value down %g ",change,value-floor(value));
761      incrementNumberTimesDown();
762      addToSumDownChange(1.0e-30+value-floor(value));
763      addToSumDownDecrease(data.intDecrease_);
764#if TYPE2==0
765      addToSumDownCost(change/(1.0e-30+(value-floor(value))));
766      setDownDynamicPseudoCost(sumDownCost()/(double) numberTimesDown());
767#elif TYPE2==1
768      addToSumDownCost(change);
769      setDownDynamicPseudoCost(sumDownCost()/sumDownChange());
770#elif TYPE2==2
771      addToSumDownCost(change*TYPERATIO+(1.0-TYPERATIO)*change/(1.0e-30+(value-floor(value))));
772      setDownDynamicPseudoCost(sumDownCost()*(TYPERATIO/sumDownChange()+(1.0-TYPERATIO)/(double) numberTimesDown()));
773#endif
774    } else {
775      //printf("(down infeasible value down %g ",change,value-floor(value));
776      incrementNumberTimesDownInfeasible();
777#if INFEAS==2
778      double distanceToCutoff=0.0;
779      double objectiveValue = model->getCurrentMinimizationObjValue();
780      distanceToCutoff =  model->getCutoff()  - originalValue;
781      if (distanceToCutoff<1.0e20) 
782        change = distanceToCutoff*2.0;
783      else 
784        change = downDynamicPseudoCost()*(value-floor(value))*10.0; 
785      change = CoinMax(1.0e-12*(1.0+fabs(originalValue)),change);
786      incrementNumberTimesDown();
787      addToSumDownChange(1.0e-30+value-floor(value));
788      addToSumDownDecrease(data.intDecrease_);
789#if TYPE2==0
790      addToSumDownCost(change/(1.0e-30+(value-floor(value))));
791      setDownDynamicPseudoCost(sumDownCost()/(double) numberTimesDown());
792#elif TYPE2==1
793      addToSumDownCost(change);
794      setDownDynamicPseudoCost(sumDownCost()/sumDownChange());
795#elif TYPE2==2
796      addToSumDownCost(change*TYPERATIO+(1.0-TYPERATIO)*change/(1.0e-30+(value-floor(value))));
797      setDownDynamicPseudoCost(sumDownCost()*(TYPERATIO/sumDownChange()+(1.0-TYPERATIO)/(double) numberTimesDown()));
798#endif
799#endif
800    }
801  } else {
802    // up
803    if (feasible) {
804      //printf("(up change %g value down %g ",change,ceil(value)-value);
805      incrementNumberTimesUp();
806      addToSumUpChange(1.0e-30+ceil(value)-value);
807      addToSumUpDecrease(data.intDecrease_);
808#if TYPE2==0
809      addToSumUpCost(change/(1.0e-30+(ceil(value)-value)));
810      setUpDynamicPseudoCost(sumUpCost()/(double) numberTimesUp());
811#elif TYPE2==1
812      addToSumUpCost(change);
813      setUpDynamicPseudoCost(sumUpCost()/sumUpChange());
814#elif TYPE2==2
815      addToSumUpCost(change*TYPERATIO+(1.0-TYPERATIO)*change/(1.0e-30+(ceil(value)-value)));
816      setUpDynamicPseudoCost(sumUpCost()*(TYPERATIO/sumUpChange()+(1.0-TYPERATIO)/(double) numberTimesUp()));
817#endif
818    } else {
819      //printf("(up infeasible value down %g ",change,ceil(value)-value);
820      incrementNumberTimesUpInfeasible();
821#if INFEAS==2
822      double distanceToCutoff=0.0;
823      double objectiveValue = model->getCurrentMinimizationObjValue();
824      distanceToCutoff =  model->getCutoff()  - originalValue;
825      if (distanceToCutoff<1.0e20) 
826        change = distanceToCutoff*2.0;
827      else 
828        change = upDynamicPseudoCost()*(ceil(value)-value)*10.0; 
829      change = CoinMax(1.0e-12*(1.0+fabs(originalValue)),change);
830      incrementNumberTimesUp();
831      addToSumUpChange(1.0e-30+ceil(value)-value);
832      addToSumUpDecrease(data.intDecrease_);
833#if TYPE2==0
834      addToSumUpCost(change/(1.0e-30+(ceil(value)-value)));
835      setUpDynamicPseudoCost(sumUpCost()/(double) numberTimesUp());
836#elif TYPE2==1
837      addToSumUpCost(change);
838      setUpDynamicPseudoCost(sumUpCost()/sumUpChange());
839#elif TYPE2==2
840      addToSumUpCost(change*TYPERATIO+(1.0-TYPERATIO)*change/(1.0e-30+(ceil(value)-value)));
841      setUpDynamicPseudoCost(sumUpCost()*(TYPERATIO/sumUpChange()+(1.0-TYPERATIO)/(double) numberTimesUp()));
842#endif
843#endif
844    }
845  }
846  //print(1,0.5);
847}
848// Pass in probing information
849void 
850CbcSimpleIntegerDynamicPseudoCost::setProbingInformation(int fixedDown, int fixedUp)
851{
852  numberTimesProbingTotal_++;
853  numberTimesDownLocalFixed_ = fixedDown;
854  numberTimesDownTotalFixed_ += fixedDown;
855  numberTimesUpLocalFixed_ = fixedUp;
856  numberTimesUpTotalFixed_ += fixedUp;
857}
858// Print
859void 
860CbcSimpleIntegerDynamicPseudoCost::print(int type,double value) const
861{
862  if (!type) {
863    double meanDown =0.0;
864    double devDown =0.0;
865    if (numberTimesDown_) {
866      meanDown = sumDownCost_/(double) numberTimesDown_;
867      devDown = meanDown*meanDown + sumDownCostSquared_ - 
868        2.0*meanDown*sumDownCost_;
869      if (devDown>=0.0)
870        devDown = sqrt(devDown);
871    }
872    double meanUp =0.0;
873    double devUp =0.0;
874    if (numberTimesUp_) {
875      meanUp = sumUpCost_/(double) numberTimesUp_;
876      devUp = meanUp*meanUp + sumUpCostSquared_ - 
877        2.0*meanUp*sumUpCost_;
878      if (devUp>=0.0)
879        devUp = sqrt(devUp);
880    }
881#if 0
882    printf("%d down %d times (%d inf) mean %g (dev %g) up %d times (%d inf) mean %g (dev %g)\n",
883           columnNumber_,
884           numberTimesDown_,numberTimesDownInfeasible_,meanDown,devDown,
885           numberTimesUp_,numberTimesUpInfeasible_,meanUp,devUp);
886#else
887    printf("%d down %d times (%d inf) mean %g  up %d times (%d inf) mean %g\n",
888           columnNumber_,
889           numberTimesDown_,numberTimesDownInfeasible_,meanDown,
890           numberTimesUp_,numberTimesUpInfeasible_,meanUp);
891#endif
892  } else {
893    const double * upper = model_->getCbcColUpper();
894    double integerTolerance = 
895      model_->getDblParam(CbcModel::CbcIntegerTolerance);
896    double below = floor(value+integerTolerance);
897    double above = below+1.0;
898    if (above>upper[columnNumber_]) {
899      above=below;
900      below = above -1;
901    }
902    double objectiveValue = model_->getCurrentMinimizationObjValue();
903    double distanceToCutoff =  model_->getCutoff() - objectiveValue;
904    if (distanceToCutoff<1.0e20) 
905      distanceToCutoff *= 10.0;
906    else 
907      distanceToCutoff = 1.0e2 + fabs(objectiveValue);
908    double sum;
909    int number;
910    double downCost = CoinMax(value-below,0.0);
911    double downCost0 = downCost*downDynamicPseudoCost_;
912    sum = sumDownCost();
913    number = numberTimesDown();
914    sum += numberTimesDownInfeasible()*(distanceToCutoff/(downCost+1.0e-12));
915    if (number>0)
916      downCost *= sum / (double) number;
917    else
918      downCost  *=  downDynamicPseudoCost_;
919    double upCost = CoinMax((above-value),0.0);
920    double upCost0 = upCost*upDynamicPseudoCost_;
921    sum = sumUpCost();
922    number = numberTimesUp();
923    sum += numberTimesUpInfeasible()*(distanceToCutoff/(upCost+1.0e-12));
924    if (number>0)
925      upCost *= sum / (double) number;
926    else
927      upCost  *=  upDynamicPseudoCost_;
928    printf("%d down %d times %g (est %g)  up %d times %g (est %g)\n",
929           columnNumber_,
930           numberTimesDown_,downCost,downCost0,
931           numberTimesUp_,upCost,upCost0);
932  }
933}
934
935// Default Constructor
936CbcDynamicPseudoCostBranchingObject::CbcDynamicPseudoCostBranchingObject()
937  :CbcIntegerBranchingObject()
938{
939  changeInGuessed_=1.0e-5;
940  object_=NULL;
941}
942
943// Useful constructor
944CbcDynamicPseudoCostBranchingObject::CbcDynamicPseudoCostBranchingObject (CbcModel * model, 
945                                                                          int variable, 
946                                                                          int way , double value,
947                                       CbcSimpleIntegerDynamicPseudoCost * object) 
948  :CbcIntegerBranchingObject(model,variable,way,value)
949{
950  changeInGuessed_=1.0e-5;
951  object_=object;
952}
953// Useful constructor for fixing
954CbcDynamicPseudoCostBranchingObject::CbcDynamicPseudoCostBranchingObject (CbcModel * model, 
955                                                      int variable, int way,
956                                                      double lowerValue, 
957                                                      double upperValue)
958  :CbcIntegerBranchingObject(model,variable,way,lowerValue)
959{
960  changeInGuessed_=1.0e100;
961  object_=NULL;
962}
963 
964
965// Copy constructor
966CbcDynamicPseudoCostBranchingObject::CbcDynamicPseudoCostBranchingObject ( 
967                                 const CbcDynamicPseudoCostBranchingObject & rhs)
968  :CbcIntegerBranchingObject(rhs)
969{
970  changeInGuessed_ = rhs.changeInGuessed_;
971  object_=rhs.object_;
972}
973
974// Assignment operator
975CbcDynamicPseudoCostBranchingObject & 
976CbcDynamicPseudoCostBranchingObject::operator=( const CbcDynamicPseudoCostBranchingObject& rhs)
977{
978  if (this != &rhs) {
979    CbcIntegerBranchingObject::operator=(rhs);
980    changeInGuessed_ = rhs.changeInGuessed_;
981    object_=rhs.object_;
982  }
983  return *this;
984}
985CbcBranchingObject * 
986CbcDynamicPseudoCostBranchingObject::clone() const
987{ 
988  return (new CbcDynamicPseudoCostBranchingObject(*this));
989}
990
991
992// Destructor
993CbcDynamicPseudoCostBranchingObject::~CbcDynamicPseudoCostBranchingObject ()
994{
995}
996
997/*
998  Perform a branch by adjusting the bounds of the specified variable. Note
999  that each arm of the branch advances the object to the next arm by
1000  advancing the value of way_.
1001
1002  Providing new values for the variable's lower and upper bounds for each
1003  branching direction gives a little bit of additional flexibility and will
1004  be easily extensible to multi-way branching.
1005  Returns change in guessed objective on next branch
1006*/
1007double
1008CbcDynamicPseudoCostBranchingObject::branch()
1009{
1010  CbcIntegerBranchingObject::branch();
1011  return changeInGuessed_;
1012}
1013/* Some branchingObjects may claim to be able to skip
1014   strong branching.  If so they have to fill in CbcStrongInfo.
1015   The object mention in incoming CbcStrongInfo must match.
1016   Returns nonzero if skip is wanted */
1017int 
1018CbcDynamicPseudoCostBranchingObject::fillStrongInfo( CbcStrongInfo & info)
1019{
1020  assert (object_);
1021  assert (info.possibleBranch==this);
1022  if (object_->numberTimesUp()<object_->numberBeforeTrust()||
1023      object_->numberTimesDown()<object_->numberBeforeTrust()) {
1024    return 0;
1025  } else {
1026    info.upMovement = object_->upDynamicPseudoCost()*(ceil(value_)-value_);
1027    info.downMovement = object_->downDynamicPseudoCost()*(value_-floor(value_));
1028    info.numIntInfeasUp  -= (int) (object_->sumUpDecrease()/
1029                                   ((double) object_->numberTimesUp()));
1030    info.numObjInfeasUp = 0;
1031    info.finishedUp = false;
1032    info.numItersUp = 0;
1033    info.numIntInfeasDown  -= (int) (object_->sumDownDecrease()/
1034                                   ((double) object_->numberTimesDown()));
1035    info.numObjInfeasDown = 0;
1036    info.finishedDown = false;
1037    info.numItersDown = 0;
1038    info.fix =0;
1039    return 1;
1040  }
1041}
1042 
1043// Default Constructor
1044CbcBranchDynamicDecision::CbcBranchDynamicDecision()
1045  :CbcBranchDecision()
1046{
1047  bestCriterion_ = 0.0;
1048  bestChangeUp_ = 0.0;
1049  bestNumberUp_ = 0;
1050  bestChangeDown_ = 0.0;
1051  bestNumberDown_ = 0;
1052  bestObject_ = NULL;
1053}
1054
1055// Copy constructor
1056CbcBranchDynamicDecision::CbcBranchDynamicDecision (
1057                                    const CbcBranchDynamicDecision & rhs)
1058  :CbcBranchDecision()
1059{
1060  bestCriterion_ = rhs.bestCriterion_;
1061  bestChangeUp_ = rhs.bestChangeUp_;
1062  bestNumberUp_ = rhs.bestNumberUp_;
1063  bestChangeDown_ = rhs.bestChangeDown_;
1064  bestNumberDown_ = rhs.bestNumberDown_;
1065  bestObject_ = rhs.bestObject_;
1066}
1067
1068CbcBranchDynamicDecision::~CbcBranchDynamicDecision()
1069{
1070}
1071
1072// Clone
1073CbcBranchDecision * 
1074CbcBranchDynamicDecision::clone() const
1075{
1076  return new CbcBranchDynamicDecision(*this);
1077}
1078
1079// Initialize i.e. before start of choosing at a node
1080void 
1081CbcBranchDynamicDecision::initialize(CbcModel * model)
1082{
1083  bestCriterion_ = 0.0;
1084  bestChangeUp_ = 0.0;
1085  bestNumberUp_ = 0;
1086  bestChangeDown_ = 0.0;
1087  bestNumberDown_ = 0;
1088  bestObject_ = NULL;
1089}
1090
1091/* Saves a clone of current branching object.  Can be used to update
1092      information on object causing branch - after branch */
1093void 
1094CbcBranchDynamicDecision::saveBranchingObject(OsiBranchingObject * object) 
1095{
1096  OsiBranchingObject * obj = object->clone();
1097  CbcDynamicPseudoCostBranchingObject * branchingObject =
1098    dynamic_cast<CbcDynamicPseudoCostBranchingObject *>(obj);
1099  assert (branchingObject);
1100  object_=branchingObject;
1101}
1102/* Pass in information on branch just done.
1103   assumes object can get information from solver */
1104void 
1105CbcBranchDynamicDecision::updateInformation(OsiSolverInterface * solver,
1106                                            const CbcNode * node)
1107{
1108  assert (object_);
1109  const CbcModel * model = object_->model();
1110  double originalValue=node->objectiveValue();
1111  int originalUnsatisfied = node->numberUnsatisfied();
1112  double objectiveValue = solver->getObjValue()*model->getObjSense();
1113  int unsatisfied=0;
1114  int i;
1115  int numberIntegers = model->numberIntegers();;
1116  const double * solution = solver->getColSolution();
1117  //const double * lower = solver->getColLower();
1118  //const double * upper = solver->getColUpper();
1119  CbcDynamicPseudoCostBranchingObject * branchingObject =
1120    dynamic_cast<CbcDynamicPseudoCostBranchingObject *>(object_);
1121  if (!branchingObject) {
1122    delete object_;
1123    object_=NULL;
1124    return;
1125  }
1126  CbcSimpleIntegerDynamicPseudoCost *  object = branchingObject->object();
1127  double change = CoinMax(0.0,objectiveValue-originalValue);
1128  // probably should also ignore if stopped
1129  int iStatus;
1130  if (solver->isProvenOptimal())
1131    iStatus=0; // optimal
1132  else if (solver->isIterationLimitReached()
1133           &&!solver->isDualObjectiveLimitReached())
1134    iStatus=2; // unknown
1135  else
1136    iStatus=1; // infeasible
1137
1138  bool feasible = iStatus!=1;
1139  if (feasible) {
1140    double integerTolerance = 
1141      model->getDblParam(CbcModel::CbcIntegerTolerance);
1142    const int * integerVariable = model->integerVariable();
1143    for (i=0;i<numberIntegers;i++) {
1144      int j=integerVariable[i];
1145      double value = solution[j];
1146      double nearest = floor(value+0.5);
1147      if (fabs(value-nearest)>integerTolerance) 
1148        unsatisfied++;
1149    }
1150  }
1151  int way = object_->way();
1152  double value = object_->value();
1153  //#define TYPE2 1
1154  //#define TYPERATIO 0.9
1155  if (way<0) {
1156    // down
1157    if (feasible) {
1158      //printf("(down change %g value down %g ",change,value-floor(value));
1159      object->incrementNumberTimesDown();
1160      object->addToSumDownChange(1.0e-30+value-floor(value));
1161      object->addToSumDownDecrease(originalUnsatisfied-unsatisfied);
1162#if TYPE2==0
1163      object->addToSumDownCost(change/(1.0e-30+(value-floor(value))));
1164      object->setDownDynamicPseudoCost(object->sumDownCost()/(double) object->numberTimesDown());
1165#elif TYPE2==1
1166      object->addToSumDownCost(change);
1167      object->setDownDynamicPseudoCost(object->sumDownCost()/object->sumDownChange());
1168#elif TYPE2==2
1169      object->addToSumDownCost(change*TYPERATIO+(1.0-TYPERATIO)*change/(1.0e-30+(value-floor(value))));
1170      object->setDownDynamicPseudoCost(object->sumDownCost()*(TYPERATIO/object->sumDownChange()+(1.0-TYPERATIO)/(double) object->numberTimesDown()));
1171#endif
1172    } else {
1173      //printf("(down infeasible value down %g ",change,value-floor(value));
1174      object->incrementNumberTimesDownInfeasible();
1175#if INFEAS==2
1176      double distanceToCutoff=0.0;
1177      double objectiveValue = model->getCurrentMinimizationObjValue();
1178      distanceToCutoff =  model->getCutoff()  - originalValue;
1179      if (distanceToCutoff<1.0e20) 
1180        change = distanceToCutoff*2.0;
1181      else 
1182        change = object->downDynamicPseudoCost()*(value-floor(value))*10.0; 
1183      change = CoinMax(1.0e-12*(1.0+fabs(originalValue)),change);
1184      object->incrementNumberTimesDown();
1185      object->addToSumDownChange(1.0e-30+value-floor(value));
1186      object->addToSumDownDecrease(originalUnsatisfied-unsatisfied);
1187#if TYPE2==0
1188      object->addToSumDownCost(change/(1.0e-30+(value-floor(value))));
1189      object->setDownDynamicPseudoCost(object->sumDownCost()/(double) object->numberTimesDown());
1190#elif TYPE2==1
1191      object->addToSumDownCost(change);
1192      object->setDownDynamicPseudoCost(object->sumDownCost()/object->sumDownChange());
1193#elif TYPE2==2
1194      object->addToSumDownCost(change*TYPERATIO+(1.0-TYPERATIO)*change/(1.0e-30+(value-floor(value))));
1195      object->setDownDynamicPseudoCost(object->sumDownCost()*(TYPERATIO/object->sumDownChange()+(1.0-TYPERATIO)/(double) object->numberTimesDown()));
1196#endif
1197#endif
1198    }
1199  } else {
1200    // up
1201    if (feasible) {
1202      //printf("(up change %g value down %g ",change,ceil(value)-value);
1203      object->incrementNumberTimesUp();
1204      object->addToSumUpChange(1.0e-30+ceil(value)-value);
1205      object->addToSumUpDecrease(unsatisfied-originalUnsatisfied);
1206#if TYPE2==0
1207      object->addToSumUpCost(change/(1.0e-30+(ceil(value)-value)));
1208      object->setUpDynamicPseudoCost(object->sumUpCost()/(double) object->numberTimesUp());
1209#elif TYPE2==1
1210      object->addToSumUpCost(change);
1211      object->setUpDynamicPseudoCost(object->sumUpCost()/object->sumUpChange());
1212#elif TYPE2==2
1213      object->addToSumUpCost(change*TYPERATIO+(1.0-TYPERATIO)*change/(1.0e-30+(ceil(value)-value)));
1214      object->setUpDynamicPseudoCost(object->sumUpCost()*(TYPERATIO/object->sumUpChange()+(1.0-TYPERATIO)/(double) object->numberTimesUp()));
1215#endif
1216    } else {
1217      //printf("(up infeasible value down %g ",change,ceil(value)-value);
1218      object->incrementNumberTimesUpInfeasible();
1219#if INFEAS==2
1220      double distanceToCutoff=0.0;
1221      double objectiveValue = model->getCurrentMinimizationObjValue();
1222      distanceToCutoff =  model->getCutoff()  - originalValue;
1223      if (distanceToCutoff<1.0e20) 
1224        change = distanceToCutoff*2.0;
1225      else 
1226        change = object->upDynamicPseudoCost()*(ceil(value)-value)*10.0; 
1227      change = CoinMax(1.0e-12*(1.0+fabs(originalValue)),change);
1228      object->incrementNumberTimesUp();
1229      object->addToSumUpChange(1.0e-30+ceil(value)-value);
1230      object->addToSumUpDecrease(unsatisfied-originalUnsatisfied);
1231#if TYPE2==0
1232      object->addToSumUpCost(change/(1.0e-30+(ceil(value)-value)));
1233      object->setUpDynamicPseudoCost(object->sumUpCost()/(double) object->numberTimesUp());
1234#elif TYPE2==1
1235      object->addToSumUpCost(change);
1236      object->setUpDynamicPseudoCost(object->sumUpCost()/object->sumUpChange());
1237#elif TYPE2==2
1238      object->addToSumUpCost(change*TYPERATIO+(1.0-TYPERATIO)*change/(1.0e-30+(ceil(value)-value)));
1239      object->setUpDynamicPseudoCost(object->sumUpCost()*(TYPERATIO/object->sumUpChange()+(1.0-TYPERATIO)/(double) object->numberTimesUp()));
1240#endif
1241#endif
1242    }
1243  }
1244  //object->print(1,0.5);
1245  delete object_;
1246  object_=NULL;
1247}
1248
1249/*
1250  Simple dynamic decision algorithm. Compare based on infeasibility (numInfUp,
1251  numInfDown) until a solution is found by search, then switch to change in
1252  objective (changeUp, changeDown). Note that bestSoFar is remembered in
1253  bestObject_, so the parameter bestSoFar is unused.
1254*/
1255
1256int
1257CbcBranchDynamicDecision::betterBranch(CbcBranchingObject * thisOne,
1258                            CbcBranchingObject * bestSoFar,
1259                            double changeUp, int numInfUp,
1260                            double changeDown, int numInfDown)
1261{
1262  CbcModel * model = thisOne->model();
1263  int stateOfSearch = model->stateOfSearch();
1264  int betterWay=0;
1265  double value=0.0;
1266  if (!bestObject_) {
1267    bestCriterion_=-1.0;
1268    bestNumberUp_=COIN_INT_MAX;
1269    bestNumberDown_=COIN_INT_MAX;
1270  }
1271  if (stateOfSearch<=1&&thisOne->model()->currentNode()->depth()>=8) {
1272#define TRY_STUFF 1
1273#ifdef TRY_STUFF
1274    // before solution - choose smallest number
1275    // could add in depth as well
1276    int bestNumber = CoinMin(bestNumberUp_,bestNumberDown_);
1277    if (numInfUp<numInfDown) {
1278      if (numInfUp<bestNumber) {
1279        betterWay = 1;
1280      } else if (numInfUp==bestNumber) {
1281        if (changeUp<bestChangeUp_)
1282          betterWay=1;
1283      }
1284    } else if (numInfUp>numInfDown) {
1285      if (numInfDown<bestNumber) {
1286        betterWay = -1;
1287      } else if (numInfDown==bestNumber) {
1288        if (changeDown<bestChangeDown_)
1289          betterWay=-1;
1290      }
1291    } else {
1292      // up and down have same number
1293      bool better=false;
1294      if (numInfUp<bestNumber) {
1295        better=true;
1296      } else if (numInfUp==bestNumber) {
1297        if (CoinMin(changeUp,changeDown)<CoinMin(bestChangeUp_,bestChangeDown_)-1.0e-5)
1298          better=true;;
1299      }
1300      if (better) {
1301        // see which way
1302        if (changeUp<=changeDown)
1303          betterWay=1;
1304        else
1305          betterWay=-1;
1306      }
1307    }
1308    if (betterWay) {
1309      value = CoinMin(numInfUp,numInfDown);
1310    }
1311#else
1312    // got a solution
1313    double minValue = CoinMin(changeDown,changeUp);
1314    double maxValue = CoinMax(changeDown,changeUp);
1315    value = WEIGHT_BEFORE*minValue + (1.0-WEIGHT_BEFORE)*maxValue;
1316    if (value>bestCriterion_+1.0e-8) {
1317      if (changeUp<=changeDown) {
1318        betterWay=1;
1319      } else {
1320        betterWay=-1;
1321      }
1322    }
1323#endif
1324  } else {
1325#if TRY_STUFF > 1
1326    // Get current number of infeasibilities, cutoff and current objective
1327    CbcNode * node = model->currentNode();
1328    int numberUnsatisfied = node->numberUnsatisfied();
1329    double cutoff = model->getCutoff();
1330    double objectiveValue = node->objectiveValue();
1331#endif
1332    // got a solution
1333    double minValue = CoinMin(changeDown,changeUp);
1334    double maxValue = CoinMax(changeDown,changeUp);
1335    // Reduce
1336#ifdef TRY_STUFF
1337    //maxValue = CoinMin(maxValue,minValue*4.0);
1338#else
1339    maxValue = CoinMin(maxValue,minValue*2.0);
1340#endif
1341    value = WEIGHT_AFTER*minValue + (1.0-WEIGHT_AFTER)*maxValue;
1342    double useValue = value;
1343    double useBest = bestCriterion_;
1344#if TRY_STUFF>1
1345    int thisNumber = CoinMin(numInfUp,numInfDown);
1346    int bestNumber = CoinMin(bestNumberUp_,bestNumberDown_);
1347    double distance = cutoff-objectiveValue;
1348    assert (distance>=0.0);
1349    if (useValue+0.1*distance>useBest&&useValue*1.1>useBest&&
1350        useBest+0.1*distance>useValue&&useBest*1.1>useValue) {
1351      // not much in it - look at unsatisfied
1352      if (thisNumber<numberUnsatisfied||bestNumber<numberUnsatisfied) {
1353        double perInteger = distance/ ((double) numberUnsatisfied);
1354        useValue += thisNumber*perInteger;
1355        useBest += bestNumber*perInteger;
1356      }
1357    }
1358#endif
1359    if (useValue>useBest+1.0e-8) {
1360      if (changeUp<=changeDown) {
1361        betterWay=1;
1362      } else {
1363        betterWay=-1;
1364      }
1365    }
1366  }
1367  if (betterWay) {
1368    // maybe change better way
1369    CbcDynamicPseudoCostBranchingObject * branchingObject =
1370      dynamic_cast<CbcDynamicPseudoCostBranchingObject *>(thisOne);
1371    if (branchingObject) {
1372      CbcSimpleIntegerDynamicPseudoCost *  object = branchingObject->object();
1373      double separator = object->upDownSeparator();
1374      if (separator>0.0) {
1375        const double * solution = thisOne->model()->testSolution();
1376        double valueVariable = solution[object->columnNumber()];
1377        betterWay = (valueVariable-floor(valueVariable)>=separator) ? 1 : -1;
1378      }
1379    }
1380    bestCriterion_ = value;
1381    bestChangeUp_ = changeUp;
1382    bestNumberUp_ = numInfUp;
1383    bestChangeDown_ = changeDown;
1384    bestNumberDown_ = numInfDown;
1385    bestObject_=thisOne;
1386    // See if user is overriding way
1387    if (thisOne->object()&&thisOne->object()->preferredWay())
1388      betterWay = thisOne->object()->preferredWay();
1389  }
1390  return betterWay;
1391}
1392/* Sets or gets best criterion so far */
1393void 
1394CbcBranchDynamicDecision::setBestCriterion(double value)
1395{ 
1396  bestCriterion_ = value;
1397}
1398double 
1399CbcBranchDynamicDecision::getBestCriterion() const
1400{ 
1401  return bestCriterion_;
1402}
Note: See TracBrowser for help on using the repository browser.