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

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

try slightly different branching

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