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

Last change on this file since 838 was 838, checked in by forrest, 14 years ago

for deterministic parallel

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 61.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#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(const 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// Updates stuff like pseudocosts before threads
380void 
381CbcSimpleIntegerDynamicPseudoCost::updateBefore(const OsiObject * rhs) 
382{
383  const CbcSimpleIntegerDynamicPseudoCost * rhsObject =
384    dynamic_cast <const CbcSimpleIntegerDynamicPseudoCost *>(rhs) ;
385  assert (rhsObject);
386  copySome(rhsObject);
387}
388  // was 1 - but that looks flakey
389#define INFEAS 1
390// Updates stuff like pseudocosts after threads finished
391void 
392CbcSimpleIntegerDynamicPseudoCost::updateAfter(const OsiObject * rhs, const OsiObject * baseObjectX) 
393{
394  const CbcSimpleIntegerDynamicPseudoCost * rhsObject =
395    dynamic_cast <const CbcSimpleIntegerDynamicPseudoCost *>(rhs) ;
396  assert (rhsObject);
397  const CbcSimpleIntegerDynamicPseudoCost * baseObject =
398    dynamic_cast <const CbcSimpleIntegerDynamicPseudoCost *>(baseObjectX) ;
399  assert (baseObject);
400  // compute current
401  double sumDown = downDynamicPseudoCost_*(numberTimesDown_+numberTimesDownInfeasible_);
402  sumDown -= baseObject->downDynamicPseudoCost_*(baseObject->numberTimesDown_+baseObject->numberTimesDownInfeasible_);
403  sumDown = CoinMax(sumDown,0.0);
404  sumDown += rhsObject->downDynamicPseudoCost_*(rhsObject->numberTimesDown_+rhsObject->numberTimesDownInfeasible_);
405  double sumUp = upDynamicPseudoCost_*(numberTimesUp_+numberTimesUpInfeasible_);
406  sumUp -= baseObject->upDynamicPseudoCost_*(baseObject->numberTimesUp_+baseObject->numberTimesUpInfeasible_);
407  sumUp += rhsObject->upDynamicPseudoCost_*(rhsObject->numberTimesUp_+rhsObject->numberTimesUpInfeasible_);
408  sumUp = CoinMax(sumUp,0.0);
409  sumDownCost_ += rhsObject->sumDownCost_-baseObject->sumDownCost_;
410  sumUpCost_ += rhsObject->sumUpCost_-baseObject->sumUpCost_;
411  sumDownChange_ += rhsObject->sumDownChange_-baseObject->sumDownChange_;
412  sumUpChange_ += rhsObject->sumUpChange_-baseObject->sumUpChange_;
413  sumDownCostSquared_ += rhsObject->sumDownCostSquared_-baseObject->sumDownCostSquared_;
414  sumUpCostSquared_ += rhsObject->sumUpCostSquared_-baseObject->sumUpCostSquared_;
415  sumDownDecrease_ += rhsObject->sumDownDecrease_-baseObject->sumDownDecrease_;
416  sumUpDecrease_ += rhsObject->sumUpDecrease_-baseObject->sumUpDecrease_;
417  lastDownCost_ += rhsObject->lastDownCost_-baseObject->lastDownCost_;
418  lastUpCost_ += rhsObject->lastUpCost_-baseObject->lastUpCost_;
419  lastDownDecrease_ += rhsObject->lastDownDecrease_-baseObject->lastDownDecrease_;
420  lastUpDecrease_ += rhsObject->lastUpDecrease_-baseObject->lastUpDecrease_;
421  numberTimesDown_ += rhsObject->numberTimesDown_-baseObject->numberTimesDown_;
422  numberTimesUp_ += rhsObject->numberTimesUp_-baseObject->numberTimesUp_;
423  numberTimesDownInfeasible_ += rhsObject->numberTimesDownInfeasible_-baseObject->numberTimesDownInfeasible_;
424  numberTimesUpInfeasible_ += rhsObject->numberTimesUpInfeasible_-baseObject->numberTimesUpInfeasible_;
425  numberTimesDownLocalFixed_ += rhsObject->numberTimesDownLocalFixed_-baseObject->numberTimesDownLocalFixed_;
426  numberTimesUpLocalFixed_ += rhsObject->numberTimesUpLocalFixed_-baseObject->numberTimesUpLocalFixed_;
427  numberTimesDownTotalFixed_ += rhsObject->numberTimesDownTotalFixed_-baseObject->numberTimesDownTotalFixed_;
428  numberTimesUpTotalFixed_ += rhsObject->numberTimesUpTotalFixed_-baseObject->numberTimesUpTotalFixed_;
429  numberTimesProbingTotal_ += rhsObject->numberTimesProbingTotal_-baseObject->numberTimesProbingTotal_;
430  if (numberTimesDown_+numberTimesDownInfeasible_>0) {
431    setDownDynamicPseudoCost(sumDown/(double) (numberTimesDown_+numberTimesDownInfeasible_));
432  }
433  if (numberTimesUp_+numberTimesUpInfeasible_>0) {
434    setUpDynamicPseudoCost(sumUp/(double) (numberTimesUp_+numberTimesUpInfeasible_));
435  }
436  //printf("XX %d down %d %d %g up %d %d %g\n",columnNumber_,numberTimesDown_,numberTimesDownInfeasible_,downDynamicPseudoCost_,
437  // numberTimesUp_,numberTimesUpInfeasible_,upDynamicPseudoCost_);
438}
439// Same - returns true if contents match(ish)
440bool 
441CbcSimpleIntegerDynamicPseudoCost::same(const CbcSimpleIntegerDynamicPseudoCost * otherObject) const
442{
443  bool okay = true;
444  if (downDynamicPseudoCost_!=otherObject->downDynamicPseudoCost_)
445    okay=false;
446  if (upDynamicPseudoCost_!=otherObject->upDynamicPseudoCost_)
447    okay=false;
448  if (sumDownCost_!= otherObject->sumDownCost_)
449    okay=false;
450  if (sumUpCost_!= otherObject->sumUpCost_)
451    okay=false;
452  if (sumDownChange_!= otherObject->sumDownChange_)
453    okay=false;
454  if (sumUpChange_!= otherObject->sumUpChange_)
455    okay=false;
456  if (sumDownCostSquared_!= otherObject->sumDownCostSquared_)
457    okay=false;
458  if (sumUpCostSquared_!= otherObject->sumUpCostSquared_)
459    okay=false;
460  if (sumDownDecrease_!= otherObject->sumDownDecrease_)
461    okay=false;
462  if (sumUpDecrease_!= otherObject->sumUpDecrease_)
463    okay=false;
464  if (lastDownCost_!= otherObject->lastDownCost_)
465    okay=false;
466  if (lastUpCost_!= otherObject->lastUpCost_)
467    okay=false;
468  if (lastDownDecrease_!= otherObject->lastDownDecrease_)
469    okay=false;
470  if (lastUpDecrease_!= otherObject->lastUpDecrease_)
471    okay=false;
472  if (numberTimesDown_!= otherObject->numberTimesDown_)
473    okay=false;
474  if (numberTimesUp_!= otherObject->numberTimesUp_)
475    okay=false;
476  if (numberTimesDownInfeasible_!= otherObject->numberTimesDownInfeasible_)
477    okay=false;
478  if (numberTimesUpInfeasible_!= otherObject->numberTimesUpInfeasible_)
479    okay=false;
480  if (numberTimesDownLocalFixed_!= otherObject->numberTimesDownLocalFixed_)
481    okay=false;
482  if (numberTimesUpLocalFixed_!= otherObject->numberTimesUpLocalFixed_)
483    okay=false;
484  if (numberTimesDownTotalFixed_!= otherObject->numberTimesDownTotalFixed_)
485    okay=false;
486  if (numberTimesUpTotalFixed_!= otherObject->numberTimesUpTotalFixed_)
487    okay=false;
488  if (numberTimesProbingTotal_!= otherObject->numberTimesProbingTotal_)
489    okay=false;
490  return okay;
491}
492// Creates a branching objecty
493CbcBranchingObject * 
494CbcSimpleIntegerDynamicPseudoCost::createBranch(int way) 
495{
496  const double * solution = model_->testSolution();
497  const double * lower = model_->getCbcColLower();
498  const double * upper = model_->getCbcColUpper();
499  double value = solution[columnNumber_];
500  value = CoinMax(value, lower[columnNumber_]);
501  value = CoinMin(value, upper[columnNumber_]);
502#ifndef NDEBUG
503  double nearest = floor(value+0.5);
504  double integerTolerance = 
505    model_->getDblParam(CbcModel::CbcIntegerTolerance);
506  assert (upper[columnNumber_]>lower[columnNumber_]);
507#endif
508  if (!model_->hotstartSolution()) {
509    assert (fabs(value-nearest)>integerTolerance);
510  } else {
511    const double * hotstartSolution = model_->hotstartSolution();
512    double targetValue = hotstartSolution[columnNumber_];
513    if (way>0)
514      value = targetValue-0.1;
515    else
516      value = targetValue+0.1;
517  }
518  CbcDynamicPseudoCostBranchingObject * newObject = 
519    new CbcDynamicPseudoCostBranchingObject(model_,columnNumber_,way,
520                                            value,this);
521  double up =  upDynamicPseudoCost_*(ceil(value)-value);
522  double down =  downDynamicPseudoCost_*(value-floor(value));
523  double changeInGuessed=up-down;
524  if (way>0)
525    changeInGuessed = - changeInGuessed;
526  changeInGuessed=CoinMax(0.0,changeInGuessed);
527  //if (way>0)
528  //changeInGuessed += 1.0e8; // bias to stay up
529  newObject->setChangeInGuessed(changeInGuessed);
530  newObject->setOriginalObject(this);
531  return newObject;
532}
533/* Create an OsiSolverBranch object
534   
535This returns NULL if branch not represented by bound changes
536*/
537OsiSolverBranch * 
538CbcSimpleIntegerDynamicPseudoCost::solverBranch() const
539{
540  OsiSolverInterface * solver = model_->solver();
541  const double * solution = model_->testSolution();
542  const double * lower = solver->getColLower();
543  const double * upper = solver->getColUpper();
544  double value = solution[columnNumber_];
545  value = CoinMax(value, lower[columnNumber_]);
546  value = CoinMin(value, upper[columnNumber_]);
547  assert (upper[columnNumber_]>lower[columnNumber_]);
548#ifndef NDEBUG
549  double nearest = floor(value+0.5);
550  double integerTolerance = 
551    model_->getDblParam(CbcModel::CbcIntegerTolerance);
552  assert (fabs(value-nearest)>integerTolerance);
553#endif
554  OsiSolverBranch * branch = new OsiSolverBranch();
555  branch->addBranch(columnNumber_,value);
556  return branch;
557}
558//#define FUNNY_BRANCHING 
559// Infeasibility - large is 0.5
560double 
561CbcSimpleIntegerDynamicPseudoCost::infeasibility(int & preferredWay) const
562{
563  const double * solution = model_->testSolution();
564  const double * lower = model_->getCbcColLower();
565  const double * upper = model_->getCbcColUpper();
566#ifdef FUNNY_BRANCHING
567  const double * dj = model_->getCbcReducedCost();
568  double djValue = dj[columnNumber_];
569  lastDownDecrease_++;
570  if (djValue>1.0e-6) {
571    // wants to go down
572    if (true||lower[columnNumber_]>originalLower_) {
573      // Lower bound active
574      lastUpDecrease_++;
575      sumDownCostSquared_ += djValue;
576    }
577  } else if (djValue<-1.0e-6) {
578    // wants to go up
579    if (true||upper[columnNumber_]<originalUpper_) {
580      // Upper bound active
581      lastUpDecrease_++;
582      sumUpCostSquared_ -= djValue;
583    }
584  }
585#endif
586  if (upper[columnNumber_]==lower[columnNumber_]) {
587    // fixed
588    preferredWay=1;
589    return 0.0;
590  }
591  assert (breakEven_>0.0&&breakEven_<1.0);
592  double value = solution[columnNumber_];
593  value = CoinMax(value, lower[columnNumber_]);
594  value = CoinMin(value, upper[columnNumber_]);
595  /*printf("%d %g %g %g %g\n",columnNumber_,value,lower[columnNumber_],
596    solution[columnNumber_],upper[columnNumber_]);*/
597  double nearest = floor(value+0.5);
598  double integerTolerance = 
599    model_->getDblParam(CbcModel::CbcIntegerTolerance);
600  double below = floor(value+integerTolerance);
601  double above = below+1.0;
602  if (above>upper[columnNumber_]) {
603    above=below;
604    below = above -1;
605  }
606#if INFEAS==1
607  double distanceToCutoff=0.0;
608  double objectiveValue = model_->getCurrentMinimizationObjValue();
609  distanceToCutoff =  model_->getCutoff()  - objectiveValue;
610  if (distanceToCutoff<1.0e20) 
611    distanceToCutoff *= 10.0;
612  else 
613    distanceToCutoff = 1.0e2 + fabs(objectiveValue);
614#endif
615  double sum;
616  double number;
617  double downCost = CoinMax(value-below,0.0);
618#if TYPE2==0
619  sum = sumDownCost_;
620  number = numberTimesDown_;
621#if INFEAS==1
622  sum += numberTimesDownInfeasible_*CoinMax(distanceToCutoff/(downCost+1.0e-12),sumDownCost_);
623  number += numberTimesDownInfeasible_;
624#endif
625#elif TYPE2==1
626  sum = sumDownCost_;
627  number = sumDownChange_;
628#if INFEAS==1
629  sum += numberTimesDownInfeasible_*CoinMax(distanceToCutoff/(downCost+1.0e-12),sumDownCost_);
630  number += numberTimesDownInfeasible_;
631#endif
632#elif TYPE2==2
633  abort();
634#if INFEAS==1
635  sum += numberTimesDownInfeasible_*(distanceToCutoff/(downCost+1.0e-12));
636  number += numberTimesDownInfeasible_;
637#endif
638#endif
639  if (number>0.0)
640    downCost *= sum / number;
641  else
642    downCost  *=  downDynamicPseudoCost_;
643  double upCost = CoinMax((above-value),0.0);
644#if TYPE2==0
645  sum = sumUpCost_;
646  number = numberTimesUp_;
647#if INFEAS==1
648  sum += numberTimesUpInfeasible_*CoinMax(distanceToCutoff/(upCost+1.0e-12),sumUpCost_);
649  number += numberTimesUpInfeasible_;
650#endif
651#elif TYPE2==1
652  sum = sumUpCost_;
653  number = sumUpChange_;
654#if INFEAS==1
655  sum += numberTimesUpInfeasible_*CoinMax(distanceToCutoff/(upCost+1.0e-12),sumUpCost_);
656  number += numberTimesUpInfeasible_;
657#endif
658#elif TYPE2==1
659  abort();
660#if INFEAS==1
661  sum += numberTimesUpInfeasible_*(distanceToCutoff/(upCost+1.0e-12));
662  number += numberTimesUpInfeasible_;
663#endif
664#endif
665  if (number>0.0)
666    upCost *= sum / number;
667  else
668    upCost  *=  upDynamicPseudoCost_;
669  if (downCost>=upCost)
670    preferredWay=1;
671  else
672    preferredWay=-1;
673  // See if up down choice set
674  if (upDownSeparator_>0.0) {
675    preferredWay = (value-below>=upDownSeparator_) ? 1 : -1;
676  }
677#ifdef FUNNY_BRANCHING
678  if (fabs(value-nearest)>integerTolerance) {
679    double ratio = (100.0+lastUpDecrease_)/(100.0+lastDownDecrease_);
680    downCost *= ratio;
681    upCost *= ratio;
682    if ((lastUpDecrease_%100)==-1) 
683      printf("col %d total %d djtimes %d down %g up %g\n",
684             columnNumber_,lastDownDecrease_,lastUpDecrease_,
685             sumDownCostSquared_,sumUpCostSquared_);
686  }
687#endif
688  if (preferredWay_)
689    preferredWay=preferredWay_;
690  // weight at 1.0 is max min
691#define WEIGHT_AFTER 0.7
692#define WEIGHT_BEFORE 0.1
693  if (fabs(value-nearest)<=integerTolerance) {
694    return 0.0;
695  } else {
696    int stateOfSearch = model_->stateOfSearch()%10;
697    double returnValue=0.0;
698    double minValue = CoinMin(downCost,upCost);
699    double maxValue = CoinMax(downCost,upCost);
700    char where;
701    // was <= 10
702    //if (stateOfSearch<=1||model_->currentNode()->depth()<=-10 /* was ||maxValue>0.2*distanceToCutoff*/) {
703    if (stateOfSearch<=2) {
704      // no branching solution
705      where='i';
706      returnValue = WEIGHT_BEFORE*minValue + (1.0-WEIGHT_BEFORE)*maxValue;
707      if (0) {
708        double sum;
709        int number;
710        double downCost2 = CoinMax(value-below,0.0);
711        sum = sumDownCost_;
712        number = numberTimesDown_;
713        if (number>0)
714          downCost2 *= sum / (double) number;
715        else
716          downCost2  *=  downDynamicPseudoCost_;
717        double upCost2 = CoinMax((above-value),0.0);
718        sum = sumUpCost_;
719        number = numberTimesUp_;
720        if (number>0)
721          upCost2 *= sum / (double) number;
722        else
723          upCost2  *=  upDynamicPseudoCost_;
724        double minValue2 = CoinMin(downCost2,upCost2);
725        double maxValue2 = CoinMax(downCost2,upCost2);
726        printf("%d value %g downC %g upC %g minV %g maxV %g downC2 %g upC2 %g minV2 %g maxV2 %g\n",
727               columnNumber_,value,downCost,upCost,minValue,maxValue,
728               downCost2,upCost2,minValue2,maxValue2);
729      }
730    } else {
731      // some solution
732      where='I';
733      returnValue = WEIGHT_AFTER*minValue + (1.0-WEIGHT_AFTER)*maxValue;
734    }
735    if (numberTimesUp_<numberBeforeTrust_||
736        numberTimesDown_<numberBeforeTrust_) {
737      //if (returnValue<1.0e10)
738      //returnValue += 1.0e12;
739      //else
740      returnValue *= 1.0e3;
741      if (!numberTimesUp_&&!numberTimesDown_)
742        returnValue *= 1.0e10;
743    }
744    //if (fabs(value-0.5)<1.0e-5) {
745    //returnValue = 3.0*returnValue + 0.2;
746    //} else if (value>0.9) {
747    //returnValue = 2.0*returnValue + 0.1;
748    //}
749    if (method_==1) {
750      // probing
751      // average
752      double up=1.0e-15;
753      double down=1.0e-15;
754      if (numberTimesProbingTotal_) {
755        up += numberTimesUpTotalFixed_/((double) numberTimesProbingTotal_);
756        down += numberTimesDownTotalFixed_/((double) numberTimesProbingTotal_);
757      }
758      returnValue = 1 + 10.0*CoinMin(numberTimesDownLocalFixed_,numberTimesUpLocalFixed_) +
759        CoinMin(down,up);
760      returnValue *= 1.0e-3;
761    }
762#ifdef COIN_DEVELOP
763    History hist;
764    hist.where_=where;
765    hist.status_=' ';
766    hist.sequence_=columnNumber_;
767    hist.numberUp_=numberTimesUp_;
768    hist.numberUpInf_=numberTimesUpInfeasible_;
769    hist.sumUp_=sumUpCost_;
770    hist.upEst_=upCost;
771    hist.numberDown_=numberTimesDown_;
772    hist.numberDownInf_=numberTimesDownInfeasible_;
773    hist.sumDown_=sumDownCost_;
774    hist.downEst_=downCost;
775    if (stateOfSearch)
776      addRecord(hist);
777#endif
778    return CoinMax(returnValue,1.0e-15);
779  }
780}
781
782double
783CbcSimpleIntegerDynamicPseudoCost::infeasibility(const OsiSolverInterface * solver, const OsiBranchingInformation * info,
784                         int & preferredWay) const
785{
786  double value = info->solution_[columnNumber_];
787  value = CoinMax(value, info->lower_[columnNumber_]);
788  value = CoinMin(value, info->upper_[columnNumber_]);
789  if (info->upper_[columnNumber_]==info->lower_[columnNumber_]) {
790    // fixed
791    preferredWay=1;
792    return 0.0;
793  }
794  assert (breakEven_>0.0&&breakEven_<1.0);
795  double nearest = floor(value+0.5);
796  double integerTolerance = info->integerTolerance_; 
797  double below = floor(value+integerTolerance);
798  double above = below+1.0;
799  if (above>info->upper_[columnNumber_]) {
800    above=below;
801    below = above -1;
802  }
803#if INFEAS==1
804  double objectiveValue = info->objectiveValue_;
805  double distanceToCutoff =  info->cutoff_  - objectiveValue;
806  if (distanceToCutoff<1.0e20) 
807    distanceToCutoff *= 10.0;
808  else 
809    distanceToCutoff = 1.0e2 + fabs(objectiveValue);
810#endif
811  double sum;
812  int number;
813  double downCost = CoinMax(value-below,0.0);
814  sum = sumDownCost_;
815  number = numberTimesDown_;
816#if INFEAS==1
817  sum += numberTimesDownInfeasible_*(distanceToCutoff/(downCost+1.0e-12));
818  number += numberTimesDownInfeasible_;
819#endif
820  if (number>0)
821    downCost *= sum / (double) number;
822  else
823    downCost  *=  downDynamicPseudoCost_;
824  double upCost = CoinMax((above-value),0.0);
825  sum = sumUpCost_;
826  number = numberTimesUp_;
827#if INFEAS==1
828  sum += numberTimesUpInfeasible_*(distanceToCutoff/(upCost+1.0e-12));
829  number += numberTimesUpInfeasible_;
830#endif
831  if (number>0)
832    upCost *= sum / (double) number;
833  else
834    upCost  *=  upDynamicPseudoCost_;
835  if (downCost>=upCost)
836    preferredWay=1;
837  else
838    preferredWay=-1;
839  // See if up down choice set
840  if (upDownSeparator_>0.0) {
841    preferredWay = (value-below>=upDownSeparator_) ? 1 : -1;
842  }
843  if (preferredWay_)
844    preferredWay=preferredWay_;
845  // weight at 1.0 is max min
846  if (fabs(value-nearest)<=integerTolerance) {
847    return 0.0;
848  } else {
849    double returnValue=0.0;
850    double minValue = CoinMin(downCost,upCost);
851    double maxValue = CoinMax(downCost,upCost);
852    if (!info->numberBranchingSolutions_||info->depth_<=10/* was ||maxValue>0.2*distanceToCutoff*/) {
853      // no solution
854      returnValue = WEIGHT_BEFORE*minValue + (1.0-WEIGHT_BEFORE)*maxValue;
855    } else {
856      // some solution
857      returnValue = WEIGHT_AFTER*minValue + (1.0-WEIGHT_AFTER)*maxValue;
858    }
859    if (numberTimesUp_<numberBeforeTrust_||
860        numberTimesDown_<numberBeforeTrust_) {
861      //if (returnValue<1.0e10)
862      //returnValue += 1.0e12;
863      //else
864      returnValue *= 1.0e3;
865      if (!numberTimesUp_&&!numberTimesDown_)
866        returnValue=1.0e50;
867    }
868    //if (fabs(value-0.5)<1.0e-5) {
869    //returnValue = 3.0*returnValue + 0.2;
870    //} else if (value>0.9) {
871    //returnValue = 2.0*returnValue + 0.1;
872    //}
873    if (method_==1) {
874      // probing
875      // average
876      double up=1.0e-15;
877      double down=1.0e-15;
878      if (numberTimesProbingTotal_) {
879        up += numberTimesUpTotalFixed_/((double) numberTimesProbingTotal_);
880        down += numberTimesDownTotalFixed_/((double) numberTimesProbingTotal_);
881      }
882      returnValue = 1 + 10.0*CoinMin(numberTimesDownLocalFixed_,numberTimesUpLocalFixed_) +
883        CoinMin(down,up);
884      returnValue *= 1.0e-3;
885    }
886    return CoinMax(returnValue,1.0e-15);
887  }
888}
889// Creates a branching object
890CbcBranchingObject * 
891CbcSimpleIntegerDynamicPseudoCost::createBranch(OsiSolverInterface * solver, const OsiBranchingInformation * info, int way) 
892{
893  double value = info->solution_[columnNumber_];
894  value = CoinMax(value, info->lower_[columnNumber_]);
895  value = CoinMin(value, info->upper_[columnNumber_]);
896  assert (info->upper_[columnNumber_]>info->lower_[columnNumber_]);
897  if (!info->hotstartSolution_) {
898#ifndef NDEBUG
899    double nearest = floor(value+0.5);
900    assert (fabs(value-nearest)>info->integerTolerance_);
901#endif
902  } else {
903    double targetValue = info->hotstartSolution_[columnNumber_];
904    if (way>0)
905      value = targetValue-0.1;
906    else
907      value = targetValue+0.1;
908  }
909  CbcDynamicPseudoCostBranchingObject * newObject = 
910    new CbcDynamicPseudoCostBranchingObject(model_,columnNumber_,way,
911                                            value,this);
912  double up =  upDynamicPseudoCost_*(ceil(value)-value);
913  double down =  downDynamicPseudoCost_*(value-floor(value));
914  double changeInGuessed=up-down;
915  if (way>0)
916    changeInGuessed = - changeInGuessed;
917  changeInGuessed=CoinMax(0.0,changeInGuessed);
918  //if (way>0)
919  //changeInGuessed += 1.0e8; // bias to stay up
920  newObject->setChangeInGuessed(changeInGuessed);
921  newObject->setOriginalObject(this);
922  return newObject;
923}
924
925// Return "up" estimate
926double 
927CbcSimpleIntegerDynamicPseudoCost::upEstimate() const
928{
929  const double * solution = model_->testSolution();
930  const double * lower = model_->getCbcColLower();
931  const double * upper = model_->getCbcColUpper();
932  double value = solution[columnNumber_];
933  value = CoinMax(value, lower[columnNumber_]);
934  value = CoinMin(value, upper[columnNumber_]);
935  if (upper[columnNumber_]==lower[columnNumber_]) {
936    // fixed
937    return 0.0;
938  }
939  double integerTolerance = 
940    model_->getDblParam(CbcModel::CbcIntegerTolerance);
941  double below = floor(value+integerTolerance);
942  double above = below+1.0;
943  if (above>upper[columnNumber_]) {
944    above=below;
945    below = above -1;
946  }
947  double upCost = CoinMax((above-value)*upDynamicPseudoCost_,0.0);
948  return upCost;
949}
950// Return "down" estimate
951double 
952CbcSimpleIntegerDynamicPseudoCost::downEstimate() const
953{
954  const double * solution = model_->testSolution();
955  const double * lower = model_->getCbcColLower();
956  const double * upper = model_->getCbcColUpper();
957  double value = solution[columnNumber_];
958  value = CoinMax(value, lower[columnNumber_]);
959  value = CoinMin(value, upper[columnNumber_]);
960  if (upper[columnNumber_]==lower[columnNumber_]) {
961    // fixed
962    return 0.0;
963  }
964  double integerTolerance = 
965    model_->getDblParam(CbcModel::CbcIntegerTolerance);
966  double below = floor(value+integerTolerance);
967  double above = below+1.0;
968  if (above>upper[columnNumber_]) {
969    above=below;
970    below = above -1;
971  }
972  double downCost = CoinMax((value-below)*downDynamicPseudoCost_,0.0);
973  return downCost;
974}
975/* Pass in information on branch just done and create CbcObjectUpdateData instance.
976   If object does not need data then backward pointer will be NULL.
977   Assumes can get information from solver */
978CbcObjectUpdateData
979CbcSimpleIntegerDynamicPseudoCost::createUpdateInformation(const OsiSolverInterface * solver, 
980                                                           const CbcNode * node,
981                                                           const CbcBranchingObject * branchingObject)
982{
983  double originalValue=node->objectiveValue();
984  int originalUnsatisfied = node->numberUnsatisfied();
985  double objectiveValue = solver->getObjValue()*solver->getObjSense();
986  int unsatisfied=0;
987  int i;
988  //might be base model - doesn't matter
989  int numberIntegers = model_->numberIntegers();;
990  const double * solution = solver->getColSolution();
991  //const double * lower = solver->getColLower();
992  //const double * upper = solver->getColUpper();
993  double change = CoinMax(0.0,objectiveValue-originalValue);
994  int iStatus;
995  if (solver->isProvenOptimal())
996    iStatus=0; // optimal
997  else if (solver->isIterationLimitReached()
998           &&!solver->isDualObjectiveLimitReached())
999    iStatus=2; // unknown
1000  else
1001    iStatus=1; // infeasible
1002
1003  bool feasible = iStatus!=1;
1004  if (feasible) {
1005    double integerTolerance = 
1006      model_->getDblParam(CbcModel::CbcIntegerTolerance);
1007    const int * integerVariable = model_->integerVariable();
1008    for (i=0;i<numberIntegers;i++) {
1009      int j=integerVariable[i];
1010      double value = solution[j];
1011      double nearest = floor(value+0.5);
1012      if (fabs(value-nearest)>integerTolerance) 
1013        unsatisfied++;
1014    }
1015  }
1016  int way = branchingObject->way();
1017  way = - way; // because after branch so moved on
1018  double value = branchingObject->value();
1019  CbcObjectUpdateData newData (this, way,
1020                               change, iStatus,
1021                               originalUnsatisfied-unsatisfied,value);
1022  newData.originalObjective_ = originalValue;
1023  // Solvers know about direction
1024  double direction = solver->getObjSense();
1025  solver->getDblParam(OsiDualObjectiveLimit,newData.cutoff_);
1026  newData.cutoff_ *= direction;
1027  return newData;
1028}
1029// Update object by CbcObjectUpdateData
1030void 
1031CbcSimpleIntegerDynamicPseudoCost::updateInformation(const CbcObjectUpdateData & data)
1032{
1033  bool feasible = data.status_!=1;
1034  int way = data.way_;
1035  double value = data.branchingValue_;
1036  double change = data.change_;
1037#define TYPERATIO 0.9
1038#define MINIMUM_MOVEMENT 0.0
1039#ifdef COIN_DEVELOP
1040  History hist;
1041  hist.where_='U'; // need to tell if hot
1042#endif
1043  double movement=0.0;
1044  if (way<0) {
1045    // down
1046    movement = value-floor(value);
1047    if (feasible) {
1048#ifdef COIN_DEVELOP
1049      hist.status_='D';
1050#endif
1051      movement = CoinMax(movement,MINIMUM_MOVEMENT);
1052      //printf("(down change %g value down %g ",change,movement);
1053      incrementNumberTimesDown();
1054      addToSumDownChange(1.0e-30+movement);
1055      addToSumDownDecrease(data.intDecrease_);
1056#if TYPE2==0
1057      addToSumDownCost(change/(1.0e-30+movement));
1058      setDownDynamicPseudoCost(sumDownCost()/(double) numberTimesDown());
1059#elif TYPE2==1
1060      addToSumDownCost(change);
1061      setDownDynamicPseudoCost(sumDownCost()/sumDownChange());
1062#elif TYPE2==2
1063      addToSumDownCost(change*TYPERATIO+(1.0-TYPERATIO)*change/(1.0e-30+movement));
1064      setDownDynamicPseudoCost(sumDownCost()*(TYPERATIO/sumDownChange()+(1.0-TYPERATIO)/(double) numberTimesDown()));
1065#endif
1066    } else {
1067#ifdef COIN_DEVELOP
1068      hist.status_='d';
1069#endif
1070      //printf("(down infeasible value down %g ",change,movement);
1071      incrementNumberTimesDownInfeasible();
1072#if INFEAS==2
1073      double distanceToCutoff=0.0;
1074      double objectiveValue = model->getCurrentMinimizationObjValue();
1075      distanceToCutoff =  model->getCutoff()  - originalValue;
1076      if (distanceToCutoff<1.0e20) 
1077        change = distanceToCutoff*2.0;
1078      else 
1079        change = downDynamicPseudoCost()*movement*10.0; 
1080      change = CoinMax(1.0e-12*(1.0+fabs(originalValue)),change);
1081      incrementNumberTimesDown();
1082      addToSumDownChange(1.0e-30+movement);
1083      addToSumDownDecrease(data.intDecrease_);
1084#if TYPE2==0
1085      addToSumDownCost(change/(1.0e-30+movement));
1086      setDownDynamicPseudoCost(sumDownCost()/(double) numberTimesDown());
1087#elif TYPE2==1
1088      addToSumDownCost(change);
1089      setDownDynamicPseudoCost(sumDownCost()/sumDownChange());
1090#elif TYPE2==2
1091      addToSumDownCost(change*TYPERATIO+(1.0-TYPERATIO)*change/(1.0e-30+movement));
1092      setDownDynamicPseudoCost(sumDownCost()*(TYPERATIO/sumDownChange()+(1.0-TYPERATIO)/(double) numberTimesDown()));
1093#endif
1094#endif
1095    }
1096#if INFEAS==1
1097    double sum = sumDownCost_;
1098    int number = numberTimesDown_;
1099    double originalValue = data.originalObjective_;
1100    assert (originalValue!=COIN_DBL_MAX);
1101    double distanceToCutoff =  data.cutoff_  - originalValue;
1102    if (distanceToCutoff>1.0e20) 
1103      distanceToCutoff=10.0+fabs(originalValue);
1104    sum += numberTimesDownInfeasible_*distanceToCutoff;
1105    number += numberTimesDownInfeasible_;
1106    setDownDynamicPseudoCost(sum/(double) number);
1107#endif
1108  } else {
1109    // up
1110    movement = ceil(value)-value;
1111    if (feasible) {
1112#ifdef COIN_DEVELOP
1113      hist.status_='U';
1114#endif
1115      movement = CoinMax(movement,MINIMUM_MOVEMENT);
1116      //printf("(up change %g value down %g ",change,movement);
1117      incrementNumberTimesUp();
1118      addToSumUpChange(1.0e-30+movement);
1119      addToSumUpDecrease(data.intDecrease_);
1120#if TYPE2==0
1121      addToSumUpCost(change/(1.0e-30+movement));
1122      setUpDynamicPseudoCost(sumUpCost()/(double) numberTimesUp());
1123#elif TYPE2==1
1124      addToSumUpCost(change);
1125      setUpDynamicPseudoCost(sumUpCost()/sumUpChange());
1126#elif TYPE2==2
1127      addToSumUpCost(change*TYPERATIO+(1.0-TYPERATIO)*change/(1.0e-30+movement));
1128      setUpDynamicPseudoCost(sumUpCost()*(TYPERATIO/sumUpChange()+(1.0-TYPERATIO)/(double) numberTimesUp()));
1129#endif
1130    } else {
1131#ifdef COIN_DEVELOP
1132      hist.status_='u';
1133#endif
1134      //printf("(up infeasible value down %g ",change,movement);
1135      incrementNumberTimesUpInfeasible();
1136#if INFEAS==2
1137      double distanceToCutoff=0.0;
1138      double objectiveValue = model->getCurrentMinimizationObjValue();
1139      distanceToCutoff =  model->getCutoff()  - originalValue;
1140      if (distanceToCutoff<1.0e20) 
1141        change = distanceToCutoff*2.0;
1142      else 
1143        change = upDynamicPseudoCost()*movement*10.0; 
1144      change = CoinMax(1.0e-12*(1.0+fabs(originalValue)),change);
1145      incrementNumberTimesUp();
1146      addToSumUpChange(1.0e-30+movement);
1147      addToSumUpDecrease(data.intDecrease_);
1148#if TYPE2==0
1149      addToSumUpCost(change/(1.0e-30+movement));
1150      setUpDynamicPseudoCost(sumUpCost()/(double) numberTimesUp());
1151#elif TYPE2==1
1152      addToSumUpCost(change);
1153      setUpDynamicPseudoCost(sumUpCost()/sumUpChange());
1154#elif TYPE2==2
1155      addToSumUpCost(change*TYPERATIO+(1.0-TYPERATIO)*change/(1.0e-30+movement));
1156      setUpDynamicPseudoCost(sumUpCost()*(TYPERATIO/sumUpChange()+(1.0-TYPERATIO)/(double) numberTimesUp()));
1157#endif
1158#endif
1159    }
1160#if INFEAS==1
1161    double sum = sumUpCost_;
1162    int number = numberTimesUp_;
1163    double originalValue = data.originalObjective_;
1164    assert (originalValue!=COIN_DBL_MAX);
1165    double distanceToCutoff =  data.cutoff_  - originalValue;
1166    if (distanceToCutoff>1.0e20) 
1167      distanceToCutoff=10.0+fabs(originalValue);
1168    sum += numberTimesUpInfeasible_*distanceToCutoff;
1169    number += numberTimesUpInfeasible_;
1170    setUpDynamicPseudoCost(sum/(double) number);
1171#endif
1172  }
1173  if (data.way_<0)
1174    assert (numberTimesDown_+numberTimesDownInfeasible_>0);
1175  else
1176    assert (numberTimesUp_+numberTimesUpInfeasible_>0);
1177  assert (downDynamicPseudoCost_>=0.0&&downDynamicPseudoCost_<1.0e100);
1178  downDynamicPseudoCost_ = CoinMax(1.0e-10,downDynamicPseudoCost_);
1179  assert (upDynamicPseudoCost_>=0.0&&upDynamicPseudoCost_<1.0e100);
1180  upDynamicPseudoCost_ = CoinMax(1.0e-10,upDynamicPseudoCost_);
1181#ifdef COIN_DEVELOP
1182  hist.sequence_=columnNumber_;
1183  hist.numberUp_=numberTimesUp_;
1184  hist.numberUpInf_=numberTimesUpInfeasible_;
1185  hist.sumUp_=sumUpCost_;
1186  hist.upEst_=change;
1187  hist.numberDown_=numberTimesDown_;
1188  hist.numberDownInf_=numberTimesDownInfeasible_;
1189  hist.sumDown_=sumDownCost_;
1190  hist.downEst_=movement;
1191  addRecord(hist);
1192#endif
1193  //print(1,0.5);
1194}
1195// Pass in probing information
1196void 
1197CbcSimpleIntegerDynamicPseudoCost::setProbingInformation(int fixedDown, int fixedUp)
1198{
1199  numberTimesProbingTotal_++;
1200  numberTimesDownLocalFixed_ = fixedDown;
1201  numberTimesDownTotalFixed_ += fixedDown;
1202  numberTimesUpLocalFixed_ = fixedUp;
1203  numberTimesUpTotalFixed_ += fixedUp;
1204}
1205// Print
1206void 
1207CbcSimpleIntegerDynamicPseudoCost::print(int type,double value) const
1208{
1209  if (!type) {
1210    double meanDown =0.0;
1211    double devDown =0.0;
1212    if (numberTimesDown_) {
1213      meanDown = sumDownCost_/(double) numberTimesDown_;
1214      devDown = meanDown*meanDown + sumDownCostSquared_ - 
1215        2.0*meanDown*sumDownCost_;
1216      if (devDown>=0.0)
1217        devDown = sqrt(devDown);
1218    }
1219    double meanUp =0.0;
1220    double devUp =0.0;
1221    if (numberTimesUp_) {
1222      meanUp = sumUpCost_/(double) numberTimesUp_;
1223      devUp = meanUp*meanUp + sumUpCostSquared_ - 
1224        2.0*meanUp*sumUpCost_;
1225      if (devUp>=0.0)
1226        devUp = sqrt(devUp);
1227    }
1228#if 0
1229    printf("%d down %d times (%d inf) mean %g (dev %g) up %d times (%d inf) mean %g (dev %g)\n",
1230           columnNumber_,
1231           numberTimesDown_,numberTimesDownInfeasible_,meanDown,devDown,
1232           numberTimesUp_,numberTimesUpInfeasible_,meanUp,devUp);
1233#else
1234    printf("%d down %d times (%d inf) mean %g  up %d times (%d inf) mean %g\n",
1235           columnNumber_,
1236           numberTimesDown_,numberTimesDownInfeasible_,meanDown,
1237           numberTimesUp_,numberTimesUpInfeasible_,meanUp);
1238#endif
1239  } else {
1240    const double * upper = model_->getCbcColUpper();
1241    double integerTolerance = 
1242      model_->getDblParam(CbcModel::CbcIntegerTolerance);
1243    double below = floor(value+integerTolerance);
1244    double above = below+1.0;
1245    if (above>upper[columnNumber_]) {
1246      above=below;
1247      below = above -1;
1248    }
1249    double objectiveValue = model_->getCurrentMinimizationObjValue();
1250    double distanceToCutoff =  model_->getCutoff() - objectiveValue;
1251    if (distanceToCutoff<1.0e20) 
1252      distanceToCutoff *= 10.0;
1253    else 
1254      distanceToCutoff = 1.0e2 + fabs(objectiveValue);
1255    double sum;
1256    int number;
1257    double downCost = CoinMax(value-below,0.0);
1258    double downCost0 = downCost*downDynamicPseudoCost_;
1259    sum = sumDownCost();
1260    number = numberTimesDown();
1261    sum += numberTimesDownInfeasible()*(distanceToCutoff/(downCost+1.0e-12));
1262    if (number>0)
1263      downCost *= sum / (double) number;
1264    else
1265      downCost  *=  downDynamicPseudoCost_;
1266    double upCost = CoinMax((above-value),0.0);
1267    double upCost0 = upCost*upDynamicPseudoCost_;
1268    sum = sumUpCost();
1269    number = numberTimesUp();
1270    sum += numberTimesUpInfeasible()*(distanceToCutoff/(upCost+1.0e-12));
1271    if (number>0)
1272      upCost *= sum / (double) number;
1273    else
1274      upCost  *=  upDynamicPseudoCost_;
1275    printf("%d down %d times %g (est %g)  up %d times %g (est %g)\n",
1276           columnNumber_,
1277           numberTimesDown_,downCost,downCost0,
1278           numberTimesUp_,upCost,upCost0);
1279  }
1280}
1281
1282// Default Constructor
1283CbcDynamicPseudoCostBranchingObject::CbcDynamicPseudoCostBranchingObject()
1284  :CbcIntegerBranchingObject()
1285{
1286  changeInGuessed_=1.0e-5;
1287  object_=NULL;
1288}
1289
1290// Useful constructor
1291CbcDynamicPseudoCostBranchingObject::CbcDynamicPseudoCostBranchingObject (CbcModel * model, 
1292                                                                          int variable, 
1293                                                                          int way , double value,
1294                                       CbcSimpleIntegerDynamicPseudoCost * object) 
1295  :CbcIntegerBranchingObject(model,variable,way,value)
1296{
1297  changeInGuessed_=1.0e-5;
1298  object_=object;
1299}
1300// Useful constructor for fixing
1301CbcDynamicPseudoCostBranchingObject::CbcDynamicPseudoCostBranchingObject (CbcModel * model, 
1302                                                      int variable, int way,
1303                                                      double lowerValue, 
1304                                                      double upperValue)
1305  :CbcIntegerBranchingObject(model,variable,way,lowerValue)
1306{
1307  changeInGuessed_=1.0e100;
1308  object_=NULL;
1309}
1310 
1311
1312// Copy constructor
1313CbcDynamicPseudoCostBranchingObject::CbcDynamicPseudoCostBranchingObject ( 
1314                                 const CbcDynamicPseudoCostBranchingObject & rhs)
1315  :CbcIntegerBranchingObject(rhs)
1316{
1317  changeInGuessed_ = rhs.changeInGuessed_;
1318  object_=rhs.object_;
1319}
1320
1321// Assignment operator
1322CbcDynamicPseudoCostBranchingObject & 
1323CbcDynamicPseudoCostBranchingObject::operator=( const CbcDynamicPseudoCostBranchingObject& rhs)
1324{
1325  if (this != &rhs) {
1326    CbcIntegerBranchingObject::operator=(rhs);
1327    changeInGuessed_ = rhs.changeInGuessed_;
1328    object_=rhs.object_;
1329  }
1330  return *this;
1331}
1332CbcBranchingObject * 
1333CbcDynamicPseudoCostBranchingObject::clone() const
1334{ 
1335  return (new CbcDynamicPseudoCostBranchingObject(*this));
1336}
1337
1338
1339// Destructor
1340CbcDynamicPseudoCostBranchingObject::~CbcDynamicPseudoCostBranchingObject ()
1341{
1342}
1343
1344/*
1345  Perform a branch by adjusting the bounds of the specified variable. Note
1346  that each arm of the branch advances the object to the next arm by
1347  advancing the value of way_.
1348
1349  Providing new values for the variable's lower and upper bounds for each
1350  branching direction gives a little bit of additional flexibility and will
1351  be easily extensible to multi-way branching.
1352  Returns change in guessed objective on next branch
1353*/
1354double
1355CbcDynamicPseudoCostBranchingObject::branch()
1356{
1357  CbcIntegerBranchingObject::branch();
1358  return changeInGuessed_;
1359}
1360/* Some branchingObjects may claim to be able to skip
1361   strong branching.  If so they have to fill in CbcStrongInfo.
1362   The object mention in incoming CbcStrongInfo must match.
1363   Returns nonzero if skip is wanted */
1364int 
1365CbcDynamicPseudoCostBranchingObject::fillStrongInfo( CbcStrongInfo & info)
1366{
1367  assert (object_);
1368  assert (info.possibleBranch==this);
1369    info.upMovement = object_->upDynamicPseudoCost()*(ceil(value_)-value_);
1370    info.downMovement = object_->downDynamicPseudoCost()*(value_-floor(value_));
1371    info.numIntInfeasUp  -= (int) (object_->sumUpDecrease()/
1372                                   (1.0e-12+(double) object_->numberTimesUp()));
1373    info.numIntInfeasUp = CoinMax(info.numIntInfeasUp,0);
1374    info.numObjInfeasUp = 0;
1375    info.finishedUp = false;
1376    info.numItersUp = 0;
1377    info.numIntInfeasDown  -= (int) (object_->sumDownDecrease()/
1378                                   (1.0e-12+(double) object_->numberTimesDown()));
1379    info.numIntInfeasDown = CoinMax(info.numIntInfeasDown,0);
1380    info.numObjInfeasDown = 0;
1381    info.finishedDown = false;
1382    info.numItersDown = 0;
1383    info.fix =0;
1384  if (object_->numberTimesUp()<object_->numberBeforeTrust()||
1385      object_->numberTimesDown()<object_->numberBeforeTrust()) {
1386    return 0;
1387  } else {
1388    return 1;
1389  }
1390}
1391 
1392// Default Constructor
1393CbcBranchDynamicDecision::CbcBranchDynamicDecision()
1394  :CbcBranchDecision()
1395{
1396  bestCriterion_ = 0.0;
1397  bestChangeUp_ = 0.0;
1398  bestNumberUp_ = 0;
1399  bestChangeDown_ = 0.0;
1400  bestNumberDown_ = 0;
1401  bestObject_ = NULL;
1402}
1403
1404// Copy constructor
1405CbcBranchDynamicDecision::CbcBranchDynamicDecision (
1406                                    const CbcBranchDynamicDecision & rhs)
1407  :CbcBranchDecision()
1408{
1409  bestCriterion_ = rhs.bestCriterion_;
1410  bestChangeUp_ = rhs.bestChangeUp_;
1411  bestNumberUp_ = rhs.bestNumberUp_;
1412  bestChangeDown_ = rhs.bestChangeDown_;
1413  bestNumberDown_ = rhs.bestNumberDown_;
1414  bestObject_ = rhs.bestObject_;
1415}
1416
1417CbcBranchDynamicDecision::~CbcBranchDynamicDecision()
1418{
1419}
1420
1421// Clone
1422CbcBranchDecision * 
1423CbcBranchDynamicDecision::clone() const
1424{
1425  return new CbcBranchDynamicDecision(*this);
1426}
1427
1428// Initialize i.e. before start of choosing at a node
1429void 
1430CbcBranchDynamicDecision::initialize(CbcModel * model)
1431{
1432  bestCriterion_ = 0.0;
1433  bestChangeUp_ = 0.0;
1434  bestNumberUp_ = 0;
1435  bestChangeDown_ = 0.0;
1436  bestNumberDown_ = 0;
1437  bestObject_ = NULL;
1438#ifdef COIN_DEVELOP
1439  History hist;
1440  hist.where_='C';
1441  hist.status_='I';
1442  hist.sequence_=55555;
1443  hist.numberUp_=0;
1444  hist.numberUpInf_=0;
1445  hist.sumUp_=0.0;
1446  hist.upEst_=0.0;
1447  hist.numberDown_=0;
1448  hist.numberDownInf_=0;
1449  hist.sumDown_=0.0;
1450  hist.downEst_=0.0;
1451  addRecord(hist);
1452#endif
1453}
1454
1455/* Saves a clone of current branching object.  Can be used to update
1456      information on object causing branch - after branch */
1457void 
1458CbcBranchDynamicDecision::saveBranchingObject(OsiBranchingObject * object) 
1459{
1460  OsiBranchingObject * obj = object->clone();
1461  CbcDynamicPseudoCostBranchingObject * branchingObject =
1462    dynamic_cast<CbcDynamicPseudoCostBranchingObject *>(obj);
1463  assert (branchingObject);
1464  object_=branchingObject;
1465}
1466/* Pass in information on branch just done.
1467   assumes object can get information from solver */
1468void 
1469CbcBranchDynamicDecision::updateInformation(OsiSolverInterface * solver,
1470                                            const CbcNode * node)
1471{
1472  assert (object_);
1473  const CbcModel * model = object_->model();
1474  double originalValue=node->objectiveValue();
1475  int originalUnsatisfied = node->numberUnsatisfied();
1476  double objectiveValue = solver->getObjValue()*model->getObjSense();
1477  int unsatisfied=0;
1478  int i;
1479  int numberIntegers = model->numberIntegers();;
1480  const double * solution = solver->getColSolution();
1481  //const double * lower = solver->getColLower();
1482  //const double * upper = solver->getColUpper();
1483  CbcDynamicPseudoCostBranchingObject * branchingObject =
1484    dynamic_cast<CbcDynamicPseudoCostBranchingObject *>(object_);
1485  if (!branchingObject) {
1486    delete object_;
1487    object_=NULL;
1488    return;
1489  }
1490  CbcSimpleIntegerDynamicPseudoCost *  object = branchingObject->object();
1491  double change = CoinMax(0.0,objectiveValue-originalValue);
1492  // probably should also ignore if stopped
1493  int iStatus;
1494  if (solver->isProvenOptimal())
1495    iStatus=0; // optimal
1496  else if (solver->isIterationLimitReached()
1497           &&!solver->isDualObjectiveLimitReached())
1498    iStatus=2; // unknown
1499  else
1500    iStatus=1; // infeasible
1501
1502  bool feasible = iStatus!=1;
1503  if (feasible) {
1504    double integerTolerance = 
1505      model->getDblParam(CbcModel::CbcIntegerTolerance);
1506    const int * integerVariable = model->integerVariable();
1507    for (i=0;i<numberIntegers;i++) {
1508      int j=integerVariable[i];
1509      double value = solution[j];
1510      double nearest = floor(value+0.5);
1511      if (fabs(value-nearest)>integerTolerance) 
1512        unsatisfied++;
1513    }
1514  }
1515  int way = object_->way();
1516  double value = object_->value();
1517  //#define TYPE2 1
1518  //#define TYPERATIO 0.9
1519  if (way<0) {
1520    // down
1521    if (feasible) {
1522      double movement = value-floor(value);
1523      movement = CoinMax(movement,MINIMUM_MOVEMENT);
1524      //printf("(down change %g value down %g ",change,movement);
1525      object->incrementNumberTimesDown();
1526      object->addToSumDownChange(1.0e-30+movement);
1527      object->addToSumDownDecrease(originalUnsatisfied-unsatisfied);
1528#if TYPE2==0
1529      object->addToSumDownCost(change/(1.0e-30+movement));
1530      object->setDownDynamicPseudoCost(object->sumDownCost()/(double) object->numberTimesDown());
1531#elif TYPE2==1
1532      object->addToSumDownCost(change);
1533      object->setDownDynamicPseudoCost(object->sumDownCost()/object->sumDownChange());
1534#elif TYPE2==2
1535      object->addToSumDownCost(change*TYPERATIO+(1.0-TYPERATIO)*change/(1.0e-30+movement));
1536      object->setDownDynamicPseudoCost(object->sumDownCost()*(TYPERATIO/object->sumDownChange()+(1.0-TYPERATIO)/(double) object->numberTimesDown()));
1537#endif
1538    } else {
1539      //printf("(down infeasible value down %g ",change,movement);
1540      object->incrementNumberTimesDownInfeasible();
1541#if INFEAS==2
1542      double distanceToCutoff=0.0;
1543      double objectiveValue = model->getCurrentMinimizationObjValue();
1544      distanceToCutoff =  model->getCutoff()  - originalValue;
1545      if (distanceToCutoff<1.0e20) 
1546        change = distanceToCutoff*2.0;
1547      else 
1548        change = object->downDynamicPseudoCost()*movement*10.0; 
1549      change = CoinMax(1.0e-12*(1.0+fabs(originalValue)),change);
1550      object->incrementNumberTimesDown();
1551      object->addToSumDownChange(1.0e-30+movement);
1552      object->addToSumDownDecrease(originalUnsatisfied-unsatisfied);
1553#if TYPE2==0
1554      object->addToSumDownCost(change/(1.0e-30+movement));
1555      object->setDownDynamicPseudoCost(object->sumDownCost()/(double) object->numberTimesDown());
1556#elif TYPE2==1
1557      object->addToSumDownCost(change);
1558      object->setDownDynamicPseudoCost(object->sumDownCost()/object->sumDownChange());
1559#elif TYPE2==2
1560      object->addToSumDownCost(change*TYPERATIO+(1.0-TYPERATIO)*change/(1.0e-30+movement));
1561      object->setDownDynamicPseudoCost(object->sumDownCost()*(TYPERATIO/object->sumDownChange()+(1.0-TYPERATIO)/(double) object->numberTimesDown()));
1562#endif
1563#endif
1564    }
1565  } else {
1566    // up
1567    if (feasible) {
1568      double movement = ceil(value)-value;
1569      movement = CoinMax(movement,MINIMUM_MOVEMENT);
1570      //printf("(up change %g value down %g ",change,movement);
1571      object->incrementNumberTimesUp();
1572      object->addToSumUpChange(1.0e-30+movement);
1573      object->addToSumUpDecrease(unsatisfied-originalUnsatisfied);
1574#if TYPE2==0
1575      object->addToSumUpCost(change/(1.0e-30+movement));
1576      object->setUpDynamicPseudoCost(object->sumUpCost()/(double) object->numberTimesUp());
1577#elif TYPE2==1
1578      object->addToSumUpCost(change);
1579      object->setUpDynamicPseudoCost(object->sumUpCost()/object->sumUpChange());
1580#elif TYPE2==2
1581      object->addToSumUpCost(change*TYPERATIO+(1.0-TYPERATIO)*change/(1.0e-30+movement));
1582      object->setUpDynamicPseudoCost(object->sumUpCost()*(TYPERATIO/object->sumUpChange()+(1.0-TYPERATIO)/(double) object->numberTimesUp()));
1583#endif
1584    } else {
1585      //printf("(up infeasible value down %g ",change,movement);
1586      object->incrementNumberTimesUpInfeasible();
1587#if INFEAS==2
1588      double distanceToCutoff=0.0;
1589      double objectiveValue = model->getCurrentMinimizationObjValue();
1590      distanceToCutoff =  model->getCutoff()  - originalValue;
1591      if (distanceToCutoff<1.0e20) 
1592        change = distanceToCutoff*2.0;
1593      else 
1594        change = object->upDynamicPseudoCost()*movement*10.0; 
1595      change = CoinMax(1.0e-12*(1.0+fabs(originalValue)),change);
1596      object->incrementNumberTimesUp();
1597      object->addToSumUpChange(1.0e-30+movement);
1598      object->addToSumUpDecrease(unsatisfied-originalUnsatisfied);
1599#if TYPE2==0
1600      object->addToSumUpCost(change/(1.0e-30+movement));
1601      object->setUpDynamicPseudoCost(object->sumUpCost()/(double) object->numberTimesUp());
1602#elif TYPE2==1
1603      object->addToSumUpCost(change);
1604      object->setUpDynamicPseudoCost(object->sumUpCost()/object->sumUpChange());
1605#elif TYPE2==2
1606      object->addToSumUpCost(change*TYPERATIO+(1.0-TYPERATIO)*change/(1.0e-30+movement));
1607      object->setUpDynamicPseudoCost(object->sumUpCost()*(TYPERATIO/object->sumUpChange()+(1.0-TYPERATIO)/(double) object->numberTimesUp()));
1608#endif
1609#endif
1610    }
1611  }
1612  //object->print(1,0.5);
1613  delete object_;
1614  object_=NULL;
1615}
1616
1617/*
1618  Simple dynamic decision algorithm. Compare based on infeasibility (numInfUp,
1619  numInfDown) until a solution is found by search, then switch to change in
1620  objective (changeUp, changeDown). Note that bestSoFar is remembered in
1621  bestObject_, so the parameter bestSoFar is unused.
1622*/
1623
1624int
1625CbcBranchDynamicDecision::betterBranch(CbcBranchingObject * thisOne,
1626                            CbcBranchingObject * bestSoFar,
1627                            double changeUp, int numInfUp,
1628                            double changeDown, int numInfDown)
1629{
1630  CbcModel * model = thisOne->model();
1631  int stateOfSearch = model->stateOfSearch()%10;
1632  int betterWay=0;
1633  double value=0.0;
1634  if (!bestObject_) {
1635    bestCriterion_=-1.0;
1636    bestNumberUp_=COIN_INT_MAX;
1637    bestNumberDown_=COIN_INT_MAX;
1638  }
1639  if (stateOfSearch<=2) {
1640    //#define TRY_STUFF 1
1641#ifdef TRY_STUFF
1642    // before solution - choose smallest number
1643    // could add in depth as well
1644    int bestNumber = CoinMin(bestNumberUp_,bestNumberDown_);
1645    if (numInfUp<numInfDown) {
1646      if (numInfUp<bestNumber) {
1647        betterWay = 1;
1648      } else if (numInfUp==bestNumber) {
1649        if (changeUp<bestChangeUp_)
1650          betterWay=1;
1651      }
1652    } else if (numInfUp>numInfDown) {
1653      if (numInfDown<bestNumber) {
1654        betterWay = -1;
1655      } else if (numInfDown==bestNumber) {
1656        if (changeDown<bestChangeDown_)
1657          betterWay=-1;
1658      }
1659    } else {
1660      // up and down have same number
1661      bool better=false;
1662      if (numInfUp<bestNumber) {
1663        better=true;
1664      } else if (numInfUp==bestNumber) {
1665        if (CoinMin(changeUp,changeDown)<CoinMin(bestChangeUp_,bestChangeDown_)-1.0e-5)
1666          better=true;;
1667      }
1668      if (better) {
1669        // see which way
1670        if (changeUp<=changeDown)
1671          betterWay=1;
1672        else
1673          betterWay=-1;
1674      }
1675    }
1676    if (betterWay) {
1677      value = CoinMin(numInfUp,numInfDown);
1678    }
1679#else
1680    // use pseudo shadow prices modified by locks
1681    // test testosi
1682#if 1
1683    double objectiveValue = model->getCurrentMinimizationObjValue();
1684    double distanceToCutoff =  model->getCutoff()  - objectiveValue;
1685    if (distanceToCutoff<1.0e20) 
1686      distanceToCutoff *= 10.0;
1687    else 
1688      distanceToCutoff = 1.0e2 + fabs(objectiveValue);
1689    double continuousObjective = model->getContinuousObjective();
1690    double distanceToCutoffC =  model->getCutoff()  - continuousObjective;
1691    if (distanceToCutoffC>1.0e20) 
1692      distanceToCutoffC = 1.0e2 + fabs(objectiveValue);
1693    int numberInfC = model->getContinuousInfeasibilities();
1694    double perInf = distanceToCutoffC/((double) numberInfC);
1695    assert (perInf>0.0);
1696    //int numberIntegers = model->numberIntegers();
1697    changeDown += perInf * numInfDown;
1698    changeUp += perInf * numInfUp;
1699#endif
1700    double minValue = CoinMin(changeDown,changeUp);
1701    double maxValue = CoinMax(changeDown,changeUp);
1702    value = WEIGHT_BEFORE*minValue + (1.0-WEIGHT_BEFORE)*maxValue;
1703    if (value>bestCriterion_+1.0e-8) {
1704      if (changeUp<=changeDown) {
1705        betterWay=1;
1706      } else {
1707        betterWay=-1;
1708      }
1709    }
1710#endif
1711  } else {
1712    //#define TRY_STUFF 2
1713#if TRY_STUFF > 1
1714    // Get current number of infeasibilities, cutoff and current objective
1715    CbcNode * node = model->currentNode();
1716    int numberUnsatisfied = node->numberUnsatisfied();
1717    double cutoff = model->getCutoff();
1718    double objectiveValue = node->objectiveValue();
1719#endif
1720    // got a solution
1721    double minValue = CoinMin(changeDown,changeUp);
1722    double maxValue = CoinMax(changeDown,changeUp);
1723    // Reduce
1724#ifdef TRY_STUFF
1725    //maxValue = CoinMin(maxValue,minValue*4.0);
1726#else
1727    //maxValue = CoinMin(maxValue,minValue*2.0);
1728#endif
1729    value = WEIGHT_AFTER*minValue + (1.0-WEIGHT_AFTER)*maxValue;
1730    double useValue = value;
1731    double useBest = bestCriterion_;
1732#if TRY_STUFF>1
1733    int thisNumber = CoinMin(numInfUp,numInfDown);
1734    int bestNumber = CoinMin(bestNumberUp_,bestNumberDown_);
1735    double distance = cutoff-objectiveValue;
1736    assert (distance>=0.0);
1737    if (useValue+0.1*distance>useBest&&useValue*1.1>useBest&&
1738        useBest+0.1*distance>useValue&&useBest*1.1>useValue) {
1739      // not much in it - look at unsatisfied
1740      if (thisNumber<numberUnsatisfied||bestNumber<numberUnsatisfied) {
1741        double perInteger = distance/ ((double) numberUnsatisfied);
1742        useValue += thisNumber*perInteger;
1743        useBest += bestNumber*perInteger;
1744      }
1745    }
1746#endif
1747    if (useValue>useBest+1.0e-8) {
1748      if (changeUp<=changeDown) {
1749        betterWay=1;
1750      } else {
1751        betterWay=-1;
1752      }
1753    }
1754  }
1755#ifdef COIN_DEVELOP
1756  History hist;
1757  {
1758    CbcDynamicPseudoCostBranchingObject * branchingObject =
1759      dynamic_cast<CbcDynamicPseudoCostBranchingObject *>(thisOne);
1760    assert (branchingObject);
1761    CbcSimpleIntegerDynamicPseudoCost *  object = branchingObject->object();
1762    assert (object);
1763    hist.where_='C';
1764    hist.status_=' ';
1765    hist.sequence_=object->columnNumber();
1766    hist.numberUp_=object->numberTimesUp();
1767    hist.numberUpInf_=numInfUp;
1768    hist.sumUp_=object->sumUpCost();
1769    hist.upEst_=changeUp;
1770    hist.numberDown_=object->numberTimesDown();
1771    hist.numberDownInf_=numInfDown;
1772    hist.sumDown_=object->sumDownCost();
1773    hist.downEst_=changeDown;
1774  }
1775#endif
1776  if (betterWay) {
1777#ifdef COIN_DEVELOP
1778    hist.status_='B';
1779#endif
1780    // maybe change better way
1781    CbcDynamicPseudoCostBranchingObject * branchingObject =
1782      dynamic_cast<CbcDynamicPseudoCostBranchingObject *>(thisOne);
1783    if (branchingObject) {
1784      CbcSimpleIntegerDynamicPseudoCost *  object = branchingObject->object();
1785      double separator = object->upDownSeparator();
1786      if (separator>0.0) {
1787        const double * solution = thisOne->model()->testSolution();
1788        double valueVariable = solution[object->columnNumber()];
1789        betterWay = (valueVariable-floor(valueVariable)>=separator) ? 1 : -1;
1790      }
1791    }
1792    bestCriterion_ = value;
1793    bestChangeUp_ = changeUp;
1794    bestNumberUp_ = numInfUp;
1795    bestChangeDown_ = changeDown;
1796    bestNumberDown_ = numInfDown;
1797    bestObject_=thisOne;
1798    // See if user is overriding way
1799    if (thisOne->object()&&thisOne->object()->preferredWay())
1800      betterWay = thisOne->object()->preferredWay();
1801  }
1802#ifdef COIN_DEVELOP
1803  addRecord(hist);
1804#endif
1805  return betterWay;
1806}
1807/* Sets or gets best criterion so far */
1808void 
1809CbcBranchDynamicDecision::setBestCriterion(double value)
1810{ 
1811  bestCriterion_ = value;
1812}
1813double 
1814CbcBranchDynamicDecision::getBestCriterion() const
1815{ 
1816  return bestCriterion_;
1817}
1818#ifdef COIN_DEVELOP
1819void printHistory(const char * file)
1820{
1821  if (!numberHistory)
1822    return;
1823  FILE * fp = fopen(file,"w");
1824  assert(fp);
1825  unsigned short numberIntegers=0;
1826  int i;
1827  for (i=0;i<numberHistory;i++) {
1828    if (history[i].where_!='C'||history[i].status_!='I') 
1829      numberIntegers = CoinMax(numberIntegers,history[i].sequence_);
1830  }
1831  numberIntegers++;
1832  for (int iC=0;iC<numberIntegers;iC++) {
1833    int n=0;
1834    for (i=0;i<numberHistory;i++) {
1835      if (history[i].sequence_==iC) {
1836        if (!n)
1837          fprintf(fp,"XXX %d\n",iC);
1838        n++;
1839        fprintf(fp,"%c%c up %8d %8d %12.5f %12.5f down  %8d %8d %12.5f %12.5f\n",
1840               history[i].where_,
1841               history[i].status_,
1842               history[i].numberUp_,
1843               history[i].numberUpInf_,
1844               history[i].sumUp_,
1845               history[i].upEst_,
1846               history[i].numberDown_,
1847               history[i].numberDownInf_,
1848               history[i].sumDown_,
1849               history[i].downEst_);
1850      }
1851    }
1852  }
1853  fclose(fp);
1854}
1855#endif
Note: See TracBrowser for help on using the repository browser.