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

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

fixes

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