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

Last change on this file since 1132 was 1132, checked in by forrest, 10 years ago

chnages to try and make faster

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