source: branches/sandbox/Cbc/src/CbcBranchDynamic.cpp @ 1286

Last change on this file since 1286 was 1286, checked in by EdwinStraver, 10 years ago

Changed formatting using AStyle -A4 -p

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