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

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

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