source: stable/2.4/Cbc/src/CbcBranchDynamic.cpp @ 1443

Last change on this file since 1443 was 1443, checked in by tkr, 10 years ago

Changing assert to warning message for unit test

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