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

Last change on this file since 892 was 887, checked in by forrest, 12 years ago

change assert to test

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