source: trunk/Cbc/src/CbcSimpleIntegerDynamicPseudoCost.cpp @ 1424

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

added Lou's comments

File size: 51.2 KB
Line 
1// Edwin 11/17/2009 - carved out of CbcBranchDynamic
2#if defined(_MSC_VER)
3// Turn off compiler warning about long names
4#  pragma warning(disable:4786)
5#endif
6#include <cassert>
7#include <cstdlib>
8#include <cmath>
9#include <cfloat>
10//#define CBC_DEBUG
11//#define TRACE_ONE 19
12#include "OsiSolverInterface.hpp"
13#include "OsiSolverBranch.hpp"
14#include "CbcModel.hpp"
15#include "CbcMessage.hpp"
16#include "CbcBranchDynamic.hpp"
17#include "CoinSort.hpp"
18#include "CoinError.hpp"
19#include "CbcSimpleIntegerDynamicPseudoCost.hpp"
20#ifdef COIN_DEVELOP
21typedef struct {
22    double sumUp_;
23    double upEst_; // or change in obj in update
24    double sumDown_;
25    double downEst_; // or movement in value in update
26    int sequence_;
27    int numberUp_;
28    int numberUpInf_;
29    int numberDown_;
30    int numberDownInf_;
31    char where_;
32    char status_;
33} History;
34History * history = NULL;
35int numberHistory = 0;
36int maxHistory = 0;
37bool getHistoryStatistics_ = true;
38static void increaseHistory()
39{
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{
50    //if (!getHistoryStatistics_)
51    return false;
52    bool fromCompare = false;
53    int i;
54    for ( i = numberHistory - 1; i >= 0; i--) {
55        if (newOne.sequence_ != history[i].sequence_)
56            continue;
57        if (newOne.where_ != history[i].where_)
58            continue;
59        if (newOne.numberUp_ != history[i].numberUp_)
60            continue;
61        if (newOne.sumUp_ != history[i].sumUp_)
62            continue;
63        if (newOne.numberUpInf_ != history[i].numberUpInf_)
64            continue;
65        if (newOne.upEst_ != history[i].upEst_)
66            continue;
67        if (newOne.numberDown_ != history[i].numberDown_)
68            continue;
69        if (newOne.sumDown_ != history[i].sumDown_)
70            continue;
71        if (newOne.numberDownInf_ != history[i].numberDownInf_)
72            continue;
73        if (newOne.downEst_ != history[i].downEst_)
74            continue;
75        // If B knock out previous B
76        if (newOne.where_ == 'C') {
77            fromCompare = true;
78            if (newOne.status_ == 'B') {
79                int j;
80                for (j = i - 1; j >= 0; j--) {
81                    if (history[j].where_ == 'C') {
82                        if (history[j].status_ == 'I') {
83                            break;
84                        } else if (history[j].status_ == 'B') {
85                            history[j].status_ = ' ';
86                            break;
87                        }
88                    }
89                }
90            }
91            break;
92        }
93    }
94    if (i == -1 || fromCompare) {
95        //add
96        increaseHistory();
97        history[numberHistory++] = newOne;
98        return true;
99    } else {
100        return false;
101    }
102}
103#endif
104
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#if TYPE2==0
186    // No
187    sumUpCost_ = 0.0;
188    sumUpChange_ = 0.0;
189    numberTimesUp_ = 0;
190    sumDownCost_ = 0.0;
191    sumDownChange_ = 0.0;
192    numberTimesDown_ = 0;
193#else
194    sumUpCost_ = 1.0 * upDynamicPseudoCost_;
195    sumUpChange_ = 1.0;
196    numberTimesUp_ = 1;
197    sumDownCost_ = 1.0 * downDynamicPseudoCost_;
198    sumDownChange_ = 1.0;
199    numberTimesDown_ = 1;
200#endif
201}
202
203/** Useful constructor
204
205  Loads dynamic upper & lower bounds for the specified variable.
206*/
207CbcSimpleIntegerDynamicPseudoCost::CbcSimpleIntegerDynamicPseudoCost (CbcModel * model,
208        int iColumn, double downDynamicPseudoCost,
209        double upDynamicPseudoCost)
210        : CbcSimpleInteger(model, iColumn),
211        upDownSeparator_(-1.0),
212        sumDownCost_(0.0),
213        sumUpCost_(0.0),
214        sumDownChange_(0.0),
215        sumUpChange_(0.0),
216        downShadowPrice_(0.0),
217        upShadowPrice_(0.0),
218        sumDownDecrease_(0.0),
219        sumUpDecrease_(0.0),
220        lastDownCost_(0.0),
221        lastUpCost_(0.0),
222        lastDownDecrease_(0),
223        lastUpDecrease_(0),
224        numberTimesDown_(0),
225        numberTimesUp_(0),
226        numberTimesDownInfeasible_(0),
227        numberTimesUpInfeasible_(0),
228        numberBeforeTrust_(0),
229        numberTimesDownLocalFixed_(0),
230        numberTimesUpLocalFixed_(0),
231        numberTimesDownTotalFixed_(0.0),
232        numberTimesUpTotalFixed_(0.0),
233        numberTimesProbingTotal_(0),
234        method_(0)
235{
236    downDynamicPseudoCost_ = downDynamicPseudoCost;
237    upDynamicPseudoCost_ = upDynamicPseudoCost;
238    breakEven_ = upDynamicPseudoCost_ / (upDynamicPseudoCost_ + downDynamicPseudoCost_);
239    // so initial will have some effect
240    sumUpCost_ = 2.0 * upDynamicPseudoCost_;
241    sumUpChange_ = 2.0;
242    numberTimesUp_ = 2;
243    sumDownCost_ = 2.0 * downDynamicPseudoCost_;
244    sumDownChange_ = 2.0;
245    numberTimesDown_ = 2;
246#if TYPE2==0
247    // No
248    sumUpCost_ = 0.0;
249    sumUpChange_ = 0.0;
250    numberTimesUp_ = 0;
251    sumDownCost_ = 0.0;
252    sumDownChange_ = 0.0;
253    numberTimesDown_ = 0;
254    sumUpCost_ = 1.0e-4 * upDynamicPseudoCost_;
255    sumDownCost_ = 1.0e-4 * downDynamicPseudoCost_;
256#else
257    sumUpCost_ = 1.0 * upDynamicPseudoCost_;
258    sumUpChange_ = 1.0;
259    numberTimesUp_ = 1;
260    sumDownCost_ = 1.0 * downDynamicPseudoCost_;
261    sumDownChange_ = 1.0;
262    numberTimesDown_ = 1;
263#endif
264}
265/** Useful constructor
266
267  Loads dynamic upper & lower bounds for the specified variable.
268*/
269CbcSimpleIntegerDynamicPseudoCost::CbcSimpleIntegerDynamicPseudoCost (CbcModel * model,
270        int /*dummy*/,
271        int iColumn, double downDynamicPseudoCost,
272        double upDynamicPseudoCost)
273{
274    CbcSimpleIntegerDynamicPseudoCost(model, iColumn, downDynamicPseudoCost, upDynamicPseudoCost);
275}
276
277// Copy constructor
278CbcSimpleIntegerDynamicPseudoCost::CbcSimpleIntegerDynamicPseudoCost ( const CbcSimpleIntegerDynamicPseudoCost & rhs)
279        : CbcSimpleInteger(rhs),
280        downDynamicPseudoCost_(rhs.downDynamicPseudoCost_),
281        upDynamicPseudoCost_(rhs.upDynamicPseudoCost_),
282        upDownSeparator_(rhs.upDownSeparator_),
283        sumDownCost_(rhs.sumDownCost_),
284        sumUpCost_(rhs.sumUpCost_),
285        sumDownChange_(rhs.sumDownChange_),
286        sumUpChange_(rhs.sumUpChange_),
287        downShadowPrice_(rhs.downShadowPrice_),
288        upShadowPrice_(rhs.upShadowPrice_),
289        sumDownDecrease_(rhs.sumDownDecrease_),
290        sumUpDecrease_(rhs.sumUpDecrease_),
291        lastDownCost_(rhs.lastDownCost_),
292        lastUpCost_(rhs.lastUpCost_),
293        lastDownDecrease_(rhs.lastDownDecrease_),
294        lastUpDecrease_(rhs.lastUpDecrease_),
295        numberTimesDown_(rhs.numberTimesDown_),
296        numberTimesUp_(rhs.numberTimesUp_),
297        numberTimesDownInfeasible_(rhs.numberTimesDownInfeasible_),
298        numberTimesUpInfeasible_(rhs.numberTimesUpInfeasible_),
299        numberBeforeTrust_(rhs.numberBeforeTrust_),
300        numberTimesDownLocalFixed_(rhs.numberTimesDownLocalFixed_),
301        numberTimesUpLocalFixed_(rhs.numberTimesUpLocalFixed_),
302        numberTimesDownTotalFixed_(rhs.numberTimesDownTotalFixed_),
303        numberTimesUpTotalFixed_(rhs.numberTimesUpTotalFixed_),
304        numberTimesProbingTotal_(rhs.numberTimesProbingTotal_),
305        method_(rhs.method_)
306
307{
308}
309
310// Clone
311CbcObject *
312CbcSimpleIntegerDynamicPseudoCost::clone() const
313{
314    return new CbcSimpleIntegerDynamicPseudoCost(*this);
315}
316
317// Assignment operator
318CbcSimpleIntegerDynamicPseudoCost &
319CbcSimpleIntegerDynamicPseudoCost::operator=( const CbcSimpleIntegerDynamicPseudoCost & rhs)
320{
321    if (this != &rhs) {
322        CbcSimpleInteger::operator=(rhs);
323        downDynamicPseudoCost_ = rhs.downDynamicPseudoCost_;
324        upDynamicPseudoCost_ = rhs.upDynamicPseudoCost_;
325        upDownSeparator_ = rhs.upDownSeparator_;
326        sumDownCost_ = rhs.sumDownCost_;
327        sumUpCost_ = rhs.sumUpCost_;
328        sumDownChange_ = rhs.sumDownChange_;
329        sumUpChange_ = rhs.sumUpChange_;
330        downShadowPrice_ = rhs.downShadowPrice_;
331        upShadowPrice_ = rhs.upShadowPrice_;
332        sumDownDecrease_ = rhs.sumDownDecrease_;
333        sumUpDecrease_ = rhs.sumUpDecrease_;
334        lastDownCost_ = rhs.lastDownCost_;
335        lastUpCost_ = rhs.lastUpCost_;
336        lastDownDecrease_ = rhs.lastDownDecrease_;
337        lastUpDecrease_ = rhs.lastUpDecrease_;
338        numberTimesDown_ = rhs.numberTimesDown_;
339        numberTimesUp_ = rhs.numberTimesUp_;
340        numberTimesDownInfeasible_ = rhs.numberTimesDownInfeasible_;
341        numberTimesUpInfeasible_ = rhs.numberTimesUpInfeasible_;
342        numberBeforeTrust_ = rhs.numberBeforeTrust_;
343        numberTimesDownLocalFixed_ = rhs.numberTimesDownLocalFixed_;
344        numberTimesUpLocalFixed_ = rhs.numberTimesUpLocalFixed_;
345        numberTimesDownTotalFixed_ = rhs.numberTimesDownTotalFixed_;
346        numberTimesUpTotalFixed_ = rhs.numberTimesUpTotalFixed_;
347        numberTimesProbingTotal_ = rhs.numberTimesProbingTotal_;
348        method_ = rhs.method_;
349    }
350    return *this;
351}
352
353// Destructor
354CbcSimpleIntegerDynamicPseudoCost::~CbcSimpleIntegerDynamicPseudoCost ()
355{
356}
357// Copy some information i.e. just variable stuff
358void
359CbcSimpleIntegerDynamicPseudoCost::copySome(const CbcSimpleIntegerDynamicPseudoCost * otherObject)
360{
361    downDynamicPseudoCost_ = otherObject->downDynamicPseudoCost_;
362    upDynamicPseudoCost_ = otherObject->upDynamicPseudoCost_;
363    sumDownCost_ = otherObject->sumDownCost_;
364    sumUpCost_ = otherObject->sumUpCost_;
365    sumDownChange_ = otherObject->sumDownChange_;
366    sumUpChange_ = otherObject->sumUpChange_;
367    downShadowPrice_ = otherObject->downShadowPrice_;
368    upShadowPrice_ = otherObject->upShadowPrice_;
369    sumDownDecrease_ = otherObject->sumDownDecrease_;
370    sumUpDecrease_ = otherObject->sumUpDecrease_;
371    lastDownCost_ = otherObject->lastDownCost_;
372    lastUpCost_ = otherObject->lastUpCost_;
373    lastDownDecrease_ = otherObject->lastDownDecrease_;
374    lastUpDecrease_ = otherObject->lastUpDecrease_;
375    numberTimesDown_ = otherObject->numberTimesDown_;
376    numberTimesUp_ = otherObject->numberTimesUp_;
377    numberTimesDownInfeasible_ = otherObject->numberTimesDownInfeasible_;
378    numberTimesUpInfeasible_ = otherObject->numberTimesUpInfeasible_;
379    numberTimesDownLocalFixed_ = otherObject->numberTimesDownLocalFixed_;
380    numberTimesUpLocalFixed_ = otherObject->numberTimesUpLocalFixed_;
381    numberTimesDownTotalFixed_ = otherObject->numberTimesDownTotalFixed_;
382    numberTimesUpTotalFixed_ = otherObject->numberTimesUpTotalFixed_;
383    numberTimesProbingTotal_ = otherObject->numberTimesProbingTotal_;
384}
385// Updates stuff like pseudocosts before threads
386void
387CbcSimpleIntegerDynamicPseudoCost::updateBefore(const OsiObject * rhs)
388{
389#ifndef NDEBUG
390    const CbcSimpleIntegerDynamicPseudoCost * rhsObject =
391        dynamic_cast <const CbcSimpleIntegerDynamicPseudoCost *>(rhs) ;
392    assert (rhsObject);
393#else
394    const CbcSimpleIntegerDynamicPseudoCost * rhsObject =
395        static_cast <const CbcSimpleIntegerDynamicPseudoCost *>(rhs) ;
396#endif
397    copySome(rhsObject);
398}
399// Updates stuff like pseudocosts after threads finished
400void
401CbcSimpleIntegerDynamicPseudoCost::updateAfter(const OsiObject * rhs, const OsiObject * baseObjectX)
402{
403#ifndef NDEBUG
404    const CbcSimpleIntegerDynamicPseudoCost * rhsObject =
405        dynamic_cast <const CbcSimpleIntegerDynamicPseudoCost *>(rhs) ;
406    assert (rhsObject);
407    const CbcSimpleIntegerDynamicPseudoCost * baseObject =
408        dynamic_cast <const CbcSimpleIntegerDynamicPseudoCost *>(baseObjectX) ;
409    assert (baseObject);
410#else
411    const CbcSimpleIntegerDynamicPseudoCost * rhsObject =
412        static_cast <const CbcSimpleIntegerDynamicPseudoCost *>(rhs) ;
413    const CbcSimpleIntegerDynamicPseudoCost * baseObject =
414        static_cast <const CbcSimpleIntegerDynamicPseudoCost *>(baseObjectX) ;
415#endif
416    // compute current
417    double sumDown = downDynamicPseudoCost_ * numberTimesDown_;
418    sumDown -= baseObject->downDynamicPseudoCost_ * baseObject->numberTimesDown_;
419    sumDown = CoinMax(sumDown, 0.0);
420    sumDown += rhsObject->downDynamicPseudoCost_ * rhsObject->numberTimesDown_;
421    assert (rhsObject->numberTimesDown_ >= baseObject->numberTimesDown_);
422    assert (rhsObject->numberTimesDownInfeasible_ >= baseObject->numberTimesDownInfeasible_);
423    assert( rhsObject->sumDownCost_ >= baseObject->sumDownCost_);
424    double sumUp = upDynamicPseudoCost_ * numberTimesUp_;
425    sumUp -= baseObject->upDynamicPseudoCost_ * baseObject->numberTimesUp_;
426    sumUp = CoinMax(sumUp, 0.0);
427    sumUp += rhsObject->upDynamicPseudoCost_ * rhsObject->numberTimesUp_;
428    assert (rhsObject->numberTimesUp_ >= baseObject->numberTimesUp_);
429    assert (rhsObject->numberTimesUpInfeasible_ >= baseObject->numberTimesUpInfeasible_);
430    assert( rhsObject->sumUpCost_ >= baseObject->sumUpCost_);
431    sumDownCost_ += rhsObject->sumDownCost_ - baseObject->sumDownCost_;
432    sumUpCost_ += rhsObject->sumUpCost_ - baseObject->sumUpCost_;
433    sumDownChange_ += rhsObject->sumDownChange_ - baseObject->sumDownChange_;
434    sumUpChange_ += rhsObject->sumUpChange_ - baseObject->sumUpChange_;
435    downShadowPrice_ = 0.0;
436    upShadowPrice_ = 0.0;
437    sumDownDecrease_ += rhsObject->sumDownDecrease_ - baseObject->sumDownDecrease_;
438    sumUpDecrease_ += rhsObject->sumUpDecrease_ - baseObject->sumUpDecrease_;
439    lastDownCost_ += rhsObject->lastDownCost_ - baseObject->lastDownCost_;
440    lastUpCost_ += rhsObject->lastUpCost_ - baseObject->lastUpCost_;
441    lastDownDecrease_ += rhsObject->lastDownDecrease_ - baseObject->lastDownDecrease_;
442    lastUpDecrease_ += rhsObject->lastUpDecrease_ - baseObject->lastUpDecrease_;
443    numberTimesDown_ += rhsObject->numberTimesDown_ - baseObject->numberTimesDown_;
444    numberTimesUp_ += rhsObject->numberTimesUp_ - baseObject->numberTimesUp_;
445    numberTimesDownInfeasible_ += rhsObject->numberTimesDownInfeasible_ - baseObject->numberTimesDownInfeasible_;
446    numberTimesUpInfeasible_ += rhsObject->numberTimesUpInfeasible_ - baseObject->numberTimesUpInfeasible_;
447    numberTimesDownLocalFixed_ += rhsObject->numberTimesDownLocalFixed_ - baseObject->numberTimesDownLocalFixed_;
448    numberTimesUpLocalFixed_ += rhsObject->numberTimesUpLocalFixed_ - baseObject->numberTimesUpLocalFixed_;
449    numberTimesDownTotalFixed_ += rhsObject->numberTimesDownTotalFixed_ - baseObject->numberTimesDownTotalFixed_;
450    numberTimesUpTotalFixed_ += rhsObject->numberTimesUpTotalFixed_ - baseObject->numberTimesUpTotalFixed_;
451    numberTimesProbingTotal_ += rhsObject->numberTimesProbingTotal_ - baseObject->numberTimesProbingTotal_;
452    if (numberTimesDown_ > 0) {
453        setDownDynamicPseudoCost(sumDown / static_cast<double> (numberTimesDown_));
454    }
455    if (numberTimesUp_ > 0) {
456        setUpDynamicPseudoCost(sumUp / static_cast<double> (numberTimesUp_));
457    }
458    //printf("XX %d down %d %d %g up %d %d %g\n",columnNumber_,numberTimesDown_,numberTimesDownInfeasible_,downDynamicPseudoCost_,
459    // numberTimesUp_,numberTimesUpInfeasible_,upDynamicPseudoCost_);
460    assert (downDynamicPseudoCost_ > 1.0e-40 && upDynamicPseudoCost_ > 1.0e-40);
461}
462// Same - returns true if contents match(ish)
463bool
464CbcSimpleIntegerDynamicPseudoCost::same(const CbcSimpleIntegerDynamicPseudoCost * otherObject) const
465{
466    bool okay = true;
467    if (downDynamicPseudoCost_ != otherObject->downDynamicPseudoCost_)
468        okay = false;
469    if (upDynamicPseudoCost_ != otherObject->upDynamicPseudoCost_)
470        okay = false;
471    if (sumDownCost_ != otherObject->sumDownCost_)
472        okay = false;
473    if (sumUpCost_ != otherObject->sumUpCost_)
474        okay = false;
475    if (sumDownChange_ != otherObject->sumDownChange_)
476        okay = false;
477    if (sumUpChange_ != otherObject->sumUpChange_)
478        okay = false;
479    if (downShadowPrice_ != otherObject->downShadowPrice_)
480        okay = false;
481    if (upShadowPrice_ != otherObject->upShadowPrice_)
482        okay = false;
483    if (sumDownDecrease_ != otherObject->sumDownDecrease_)
484        okay = false;
485    if (sumUpDecrease_ != otherObject->sumUpDecrease_)
486        okay = false;
487    if (lastDownCost_ != otherObject->lastDownCost_)
488        okay = false;
489    if (lastUpCost_ != otherObject->lastUpCost_)
490        okay = false;
491    if (lastDownDecrease_ != otherObject->lastDownDecrease_)
492        okay = false;
493    if (lastUpDecrease_ != otherObject->lastUpDecrease_)
494        okay = false;
495    if (numberTimesDown_ != otherObject->numberTimesDown_)
496        okay = false;
497    if (numberTimesUp_ != otherObject->numberTimesUp_)
498        okay = false;
499    if (numberTimesDownInfeasible_ != otherObject->numberTimesDownInfeasible_)
500        okay = false;
501    if (numberTimesUpInfeasible_ != otherObject->numberTimesUpInfeasible_)
502        okay = false;
503    if (numberTimesDownLocalFixed_ != otherObject->numberTimesDownLocalFixed_)
504        okay = false;
505    if (numberTimesUpLocalFixed_ != otherObject->numberTimesUpLocalFixed_)
506        okay = false;
507    if (numberTimesDownTotalFixed_ != otherObject->numberTimesDownTotalFixed_)
508        okay = false;
509    if (numberTimesUpTotalFixed_ != otherObject->numberTimesUpTotalFixed_)
510        okay = false;
511    if (numberTimesProbingTotal_ != otherObject->numberTimesProbingTotal_)
512        okay = false;
513    return okay;
514}
515/* Create an OsiSolverBranch object
516
517This returns NULL if branch not represented by bound changes
518*/
519OsiSolverBranch *
520CbcSimpleIntegerDynamicPseudoCost::solverBranch() const
521{
522    OsiSolverInterface * solver = model_->solver();
523    const double * solution = model_->testSolution();
524    const double * lower = solver->getColLower();
525    const double * upper = solver->getColUpper();
526    double value = solution[columnNumber_];
527    value = CoinMax(value, lower[columnNumber_]);
528    value = CoinMin(value, upper[columnNumber_]);
529    assert (upper[columnNumber_] > lower[columnNumber_]);
530#ifndef NDEBUG
531    double nearest = floor(value + 0.5);
532    double integerTolerance =
533        model_->getDblParam(CbcModel::CbcIntegerTolerance);
534    assert (fabs(value - nearest) > integerTolerance);
535#endif
536    OsiSolverBranch * branch = new OsiSolverBranch();
537    branch->addBranch(columnNumber_, value);
538    return branch;
539}
540//#define FUNNY_BRANCHING
541double
542CbcSimpleIntegerDynamicPseudoCost::infeasibility(const OsiBranchingInformation * info,
543        int &preferredWay) const
544{
545    assert (downDynamicPseudoCost_ > 1.0e-40 && upDynamicPseudoCost_ > 1.0e-40);
546    const double * solution = model_->testSolution();
547    const double * lower = model_->getCbcColLower();
548    const double * upper = model_->getCbcColUpper();
549#ifdef FUNNY_BRANCHING
550    const double * dj = model_->getCbcReducedCost();
551    double djValue = dj[columnNumber_];
552    lastDownDecrease_++;
553    if (djValue > 1.0e-6) {
554        // wants to go down
555        if (true || lower[columnNumber_] > originalLower_) {
556            // Lower bound active
557            lastUpDecrease_++;
558        }
559    } else if (djValue < -1.0e-6) {
560        // wants to go up
561        if (true || upper[columnNumber_] < originalUpper_) {
562            // Upper bound active
563            lastUpDecrease_++;
564        }
565    }
566#endif
567    if (upper[columnNumber_] == lower[columnNumber_]) {
568        // fixed
569        preferredWay = 1;
570        return 0.0;
571    }
572    assert (breakEven_ > 0.0 && breakEven_ < 1.0);
573        /*
574          Find nearest integer, and integers above and below current value.
575
576          Given that we've already forced value within bounds, if
577          (current value)+(integer tolerance) > (upper bound)
578          shouldn't we declare this variable integer?
579        */
580
581    double value = solution[columnNumber_];
582    value = CoinMax(value, lower[columnNumber_]);
583    value = CoinMin(value, upper[columnNumber_]);
584    /*printf("%d %g %g %g %g\n",columnNumber_,value,lower[columnNumber_],
585      solution[columnNumber_],upper[columnNumber_]);*/
586    double nearest = floor(value + 0.5);
587    double integerTolerance =
588        model_->getDblParam(CbcModel::CbcIntegerTolerance);
589    double below = floor(value + integerTolerance);
590    double above = below + 1.0;
591    if (above > upper[columnNumber_]) {
592        above = below;
593        below = above - 1;
594    }
595#if INFEAS==1
596/*
597  Why do we inflate the distance to the cutoff by a factor of 10 for
598  values that could be considered reachable? Why do we add 100 for values
599  larger than 1e20?
600*/
601    double distanceToCutoff = 0.0;
602    double objectiveValue = model_->getCurrentMinimizationObjValue();
603    distanceToCutoff =  model_->getCutoff()  - objectiveValue;
604    if (distanceToCutoff < 1.0e20)
605        distanceToCutoff *= 10.0;
606    else
607        distanceToCutoff = 1.0e2 + fabs(objectiveValue);
608    distanceToCutoff = CoinMax(distanceToCutoff, 1.0e-12 * (1.0 + fabs(objectiveValue)));
609#endif
610    double sum;
611    double number;
612    double downCost = CoinMax(value - below, 0.0);
613#if TYPE2==0
614    sum = sumDownCost_;
615    number = numberTimesDown_;
616#if INFEAS==1
617    sum += numberTimesDownInfeasible_ * CoinMax(distanceToCutoff / (downCost + 1.0e-12), sumDownCost_);
618#endif
619#elif TYPE2==1
620    sum = sumDownCost_;
621    number = sumDownChange_;
622#if INFEAS==1
623    sum += numberTimesDownInfeasible_ * CoinMax(distanceToCutoff / (downCost + 1.0e-12), sumDownCost_);
624#endif
625#elif TYPE2==2
626    abort();
627#if INFEAS==1
628    sum += numberTimesDownInfeasible_ * (distanceToCutoff / (downCost + 1.0e-12));
629#endif
630#endif
631#if MOD_SHADOW>0
632    if (!downShadowPrice_) {
633        if (number > 0.0)
634            downCost *= sum / number;
635        else
636            downCost  *=  downDynamicPseudoCost_;
637    } else if (downShadowPrice_ > 0.0) {
638        downCost *= downShadowPrice_;
639    } else {
640        downCost *= (downDynamicPseudoCost_ - downShadowPrice_);
641    }
642#else
643    if (downShadowPrice_ <= 0.0) {
644        if (number > 0.0)
645            downCost *= sum / number;
646        else
647            downCost  *=  downDynamicPseudoCost_;
648    } else {
649        downCost *= downShadowPrice_;
650    }
651#endif
652    double upCost = CoinMax((above - value), 0.0);
653#if TYPE2==0
654    sum = sumUpCost_;
655    number = numberTimesUp_;
656#if INFEAS==1
657    sum += numberTimesUpInfeasible_ * CoinMax(distanceToCutoff / (upCost + 1.0e-12), sumUpCost_);
658#endif
659#elif TYPE2==1
660    sum = sumUpCost_;
661    number = sumUpChange_;
662#if INFEAS==1
663    sum += numberTimesUpInfeasible_ * CoinMax(distanceToCutoff / (upCost + 1.0e-12), sumUpCost_);
664#endif
665#elif TYPE2==1
666    abort();
667#if INFEAS==1
668    sum += numberTimesUpInfeasible_ * (distanceToCutoff / (upCost + 1.0e-12));
669#endif
670#endif
671#if MOD_SHADOW>0
672    if (!upShadowPrice_) {
673        if (number > 0.0)
674            upCost *= sum / number;
675        else
676            upCost  *=  upDynamicPseudoCost_;
677    } else if (upShadowPrice_ > 0.0) {
678        upCost *= upShadowPrice_;
679    } else {
680        upCost *= (upDynamicPseudoCost_ - upShadowPrice_);
681    }
682#else
683    if (upShadowPrice_ <= 0.0) {
684        if (number > 0.0)
685            upCost *= sum / number;
686        else
687            upCost  *=  upDynamicPseudoCost_;
688    } else {
689        upCost *= upShadowPrice_;
690    }
691#endif
692    if (downCost >= upCost)
693        preferredWay = 1;
694    else
695        preferredWay = -1;
696    // See if up down choice set
697    if (upDownSeparator_ > 0.0) {
698        preferredWay = (value - below >= upDownSeparator_) ? 1 : -1;
699    }
700#ifdef FUNNY_BRANCHING
701    if (fabs(value - nearest) > integerTolerance) {
702        double ratio = (100.0 + lastUpDecrease_) / (100.0 + lastDownDecrease_);
703        downCost *= ratio;
704        upCost *= ratio;
705        if ((lastUpDecrease_ % 100) == -1)
706            printf("col %d total %d djtimes %d\n",
707                   columnNumber_, lastDownDecrease_, lastUpDecrease_);
708    }
709#endif
710    if (preferredWay_)
711        preferredWay = preferredWay_;
712    if (info->hotstartSolution_) {
713        double targetValue = info->hotstartSolution_[columnNumber_];
714        if (value > targetValue)
715            preferredWay = -1;
716        else
717            preferredWay = 1;
718    }
719    if (fabs(value - nearest) <= integerTolerance) {
720        if (priority_ != -999)
721            return 0.0;
722        else
723            return 1.0e-13;
724    } else {
725        int stateOfSearch = model_->stateOfSearch() % 10;
726        double returnValue = 0.0;
727        double minValue = CoinMin(downCost, upCost);
728        double maxValue = CoinMax(downCost, upCost);
729        char where;
730        // was <= 10
731        //if (stateOfSearch<=1||model_->currentNode()->depth()<=-10 /* was ||maxValue>0.2*distanceToCutoff*/) {
732        if (stateOfSearch <= 2) {
733            // no branching solution
734            where = 'i';
735            returnValue = WEIGHT_BEFORE * minValue + (1.0 - WEIGHT_BEFORE) * maxValue;
736            if (0) {
737                double sum;
738                int number;
739                double downCost2 = CoinMax(value - below, 0.0);
740                sum = sumDownCost_;
741                number = numberTimesDown_;
742                if (number > 0)
743                    downCost2 *= sum / static_cast<double> (number);
744                else
745                    downCost2  *=  downDynamicPseudoCost_;
746                double upCost2 = CoinMax((above - value), 0.0);
747                sum = sumUpCost_;
748                number = numberTimesUp_;
749                if (number > 0)
750                    upCost2 *= sum / static_cast<double> (number);
751                else
752                    upCost2  *=  upDynamicPseudoCost_;
753                double minValue2 = CoinMin(downCost2, upCost2);
754                double maxValue2 = CoinMax(downCost2, upCost2);
755                printf("%d value %g downC %g upC %g minV %g maxV %g downC2 %g upC2 %g minV2 %g maxV2 %g\n",
756                       columnNumber_, value, downCost, upCost, minValue, maxValue,
757                       downCost2, upCost2, minValue2, maxValue2);
758            }
759        } else {
760            // some solution
761            where = 'I';
762#ifndef WEIGHT_PRODUCT
763            returnValue = WEIGHT_AFTER * minValue + (1.0 - WEIGHT_AFTER) * maxValue;
764#else
765            double minProductWeight = model_->getDblParam(CbcModel::CbcSmallChange);
766            returnValue = CoinMax(minValue, minProductWeight) * CoinMax(maxValue, minProductWeight);
767            //returnValue += minProductWeight*minValue;
768#endif
769        }
770        if (numberTimesUp_ < numberBeforeTrust_ ||
771                numberTimesDown_ < numberBeforeTrust_) {
772            //if (returnValue<1.0e10)
773            //returnValue += 1.0e12;
774            //else
775            returnValue *= 1.0e3;
776            if (!numberTimesUp_ && !numberTimesDown_)
777                returnValue *= 1.0e10;
778        }
779        //if (fabs(value-0.5)<1.0e-5) {
780        //returnValue = 3.0*returnValue + 0.2;
781        //} else if (value>0.9) {
782        //returnValue = 2.0*returnValue + 0.1;
783        //}
784        if (method_ == 1) {
785            // probing
786            // average
787            double up = 1.0e-15;
788            double down = 1.0e-15;
789            if (numberTimesProbingTotal_) {
790                up += numberTimesUpTotalFixed_ / static_cast<double> (numberTimesProbingTotal_);
791                down += numberTimesDownTotalFixed_ / static_cast<double> (numberTimesProbingTotal_);
792            }
793            returnValue = 1 + 10.0 * CoinMin(numberTimesDownLocalFixed_, numberTimesUpLocalFixed_) +
794                          CoinMin(down, up);
795            returnValue *= 1.0e-3;
796        }
797#ifdef COIN_DEVELOP
798        History hist;
799        hist.where_ = where;
800        hist.status_ = ' ';
801        hist.sequence_ = columnNumber_;
802        hist.numberUp_ = numberTimesUp_;
803        hist.numberUpInf_ = numberTimesUpInfeasible_;
804        hist.sumUp_ = sumUpCost_;
805        hist.upEst_ = upCost;
806        hist.numberDown_ = numberTimesDown_;
807        hist.numberDownInf_ = numberTimesDownInfeasible_;
808        hist.sumDown_ = sumDownCost_;
809        hist.downEst_ = downCost;
810        if (stateOfSearch)
811            addRecord(hist);
812#endif
813        return CoinMax(returnValue, 1.0e-15);
814    }
815}
816// Creates a branching object
817CbcBranchingObject *
818CbcSimpleIntegerDynamicPseudoCost::createCbcBranch(OsiSolverInterface * /*solver*/,
819        const OsiBranchingInformation * info, int way)
820{
821    double value = info->solution_[columnNumber_];
822    value = CoinMax(value, info->lower_[columnNumber_]);
823    value = CoinMin(value, info->upper_[columnNumber_]);
824    assert (info->upper_[columnNumber_] > info->lower_[columnNumber_]);
825    if (!info->hotstartSolution_ && priority_ != -999) {
826#ifndef NDEBUG
827        double nearest = floor(value + 0.5);
828        assert (fabs(value - nearest) > info->integerTolerance_);
829#endif
830    } else if (info->hotstartSolution_) {
831        double targetValue = info->hotstartSolution_[columnNumber_];
832        if (way > 0)
833            value = targetValue - 0.1;
834        else
835            value = targetValue + 0.1;
836    } else {
837        if (value <= info->lower_[columnNumber_])
838            value += 0.1;
839        else if (value >= info->upper_[columnNumber_])
840            value -= 0.1;
841    }
842    assert (value >= info->lower_[columnNumber_] &&
843            value <= info->upper_[columnNumber_]);
844    CbcDynamicPseudoCostBranchingObject * newObject =
845        new CbcDynamicPseudoCostBranchingObject(model_, columnNumber_, way,
846                                                value, this);
847    double up =  upDynamicPseudoCost_ * (ceil(value) - value);
848    double down =  downDynamicPseudoCost_ * (value - floor(value));
849    double changeInGuessed = up - down;
850    if (way > 0)
851        changeInGuessed = - changeInGuessed;
852    changeInGuessed = CoinMax(0.0, changeInGuessed);
853    //if (way>0)
854    //changeInGuessed += 1.0e8; // bias to stay up
855    newObject->setChangeInGuessed(changeInGuessed);
856    newObject->setOriginalObject(this);
857    return newObject;
858}
859
860// Return "up" estimate
861double
862CbcSimpleIntegerDynamicPseudoCost::upEstimate() const
863{
864    const double * solution = model_->testSolution();
865    const double * lower = model_->getCbcColLower();
866    const double * upper = model_->getCbcColUpper();
867    double value = solution[columnNumber_];
868    value = CoinMax(value, lower[columnNumber_]);
869    value = CoinMin(value, upper[columnNumber_]);
870    if (upper[columnNumber_] == lower[columnNumber_]) {
871        // fixed
872        return 0.0;
873    }
874    double integerTolerance =
875        model_->getDblParam(CbcModel::CbcIntegerTolerance);
876    double below = floor(value + integerTolerance);
877    double above = below + 1.0;
878    if (above > upper[columnNumber_]) {
879        above = below;
880        below = above - 1;
881    }
882    double upCost = CoinMax((above - value) * upDynamicPseudoCost_, 0.0);
883    return upCost;
884}
885// Return "down" estimate
886double
887CbcSimpleIntegerDynamicPseudoCost::downEstimate() const
888{
889    const double * solution = model_->testSolution();
890    const double * lower = model_->getCbcColLower();
891    const double * upper = model_->getCbcColUpper();
892    double value = solution[columnNumber_];
893    value = CoinMax(value, lower[columnNumber_]);
894    value = CoinMin(value, upper[columnNumber_]);
895    if (upper[columnNumber_] == lower[columnNumber_]) {
896        // fixed
897        return 0.0;
898    }
899    double integerTolerance =
900        model_->getDblParam(CbcModel::CbcIntegerTolerance);
901    double below = floor(value + integerTolerance);
902    double above = below + 1.0;
903    if (above > upper[columnNumber_]) {
904        above = below;
905        below = above - 1;
906    }
907    double downCost = CoinMax((value - below) * downDynamicPseudoCost_, 0.0);
908    return downCost;
909}
910// Set down pseudo cost
911void
912CbcSimpleIntegerDynamicPseudoCost::setDownDynamicPseudoCost(double value)
913{
914#ifdef TRACE_ONE
915    double oldDown = sumDownCost_;
916#endif
917    downDynamicPseudoCost_ = value;
918    sumDownCost_ = CoinMax(sumDownCost_, value * numberTimesDown_);
919#ifdef TRACE_ONE
920    if (columnNumber_ == TRACE_ONE) {
921        double down = downDynamicPseudoCost_ * numberTimesDown_;
922        printf("For %d sumDown %g (%d), inf (%d) - pseudo %g - sumDown was %g -> %g\n",
923               TRACE_ONE, down, numberTimesDown_,
924               numberTimesDownInfeasible_, downDynamicPseudoCost_,
925               oldDown, sumDownCost_);
926    }
927#endif
928}
929// Modify down pseudo cost in a slightly different way
930void
931CbcSimpleIntegerDynamicPseudoCost::updateDownDynamicPseudoCost(double value)
932{
933    sumDownCost_ += value;
934    numberTimesDown_++;
935    downDynamicPseudoCost_ = sumDownCost_ / static_cast<double>(numberTimesDown_);
936}
937// Set up pseudo cost
938void
939CbcSimpleIntegerDynamicPseudoCost::setUpDynamicPseudoCost(double value)
940{
941#ifdef TRACE_ONE
942    double oldUp = sumUpCost_;
943#endif
944    upDynamicPseudoCost_ = value;
945    sumUpCost_ = CoinMax(sumUpCost_, value * numberTimesUp_);
946#ifdef TRACE_ONE
947    if (columnNumber_ == TRACE_ONE) {
948        double up = upDynamicPseudoCost_ * numberTimesUp_;
949        printf("For %d sumUp %g (%d), inf (%d) - pseudo %g - sumUp was %g -> %g\n",
950               TRACE_ONE, up, numberTimesUp_,
951               numberTimesUpInfeasible_, upDynamicPseudoCost_,
952               oldUp, sumUpCost_);
953    }
954#endif
955}
956// Modify up pseudo cost in a slightly different way
957void
958CbcSimpleIntegerDynamicPseudoCost::updateUpDynamicPseudoCost(double value)
959{
960    sumUpCost_ += value;
961    numberTimesUp_++;
962    upDynamicPseudoCost_ = sumUpCost_ / static_cast<double>(numberTimesUp_);
963}
964/* Pass in information on branch just done and create CbcObjectUpdateData instance.
965   If object does not need data then backward pointer will be NULL.
966   Assumes can get information from solver */
967CbcObjectUpdateData
968CbcSimpleIntegerDynamicPseudoCost::createUpdateInformation(const OsiSolverInterface * solver,
969        const CbcNode * node,
970        const CbcBranchingObject * branchingObject)
971{
972    double originalValue = node->objectiveValue();
973    int originalUnsatisfied = node->numberUnsatisfied();
974    double objectiveValue = solver->getObjValue() * solver->getObjSense();
975    int unsatisfied = 0;
976    int i;
977    //might be base model - doesn't matter
978    int numberIntegers = model_->numberIntegers();;
979    const double * solution = solver->getColSolution();
980    //const double * lower = solver->getColLower();
981    //const double * upper = solver->getColUpper();
982    double change = CoinMax(0.0, objectiveValue - originalValue);
983    int iStatus;
984    if (solver->isProvenOptimal())
985        iStatus = 0; // optimal
986    else if (solver->isIterationLimitReached()
987             && !solver->isDualObjectiveLimitReached())
988        iStatus = 2; // unknown
989    else
990        iStatus = 1; // infeasible
991
992    bool feasible = iStatus != 1;
993    if (feasible) {
994        double integerTolerance =
995            model_->getDblParam(CbcModel::CbcIntegerTolerance);
996        const int * integerVariable = model_->integerVariable();
997        for (i = 0; i < numberIntegers; i++) {
998            int j = integerVariable[i];
999            double value = solution[j];
1000            double nearest = floor(value + 0.5);
1001            if (fabs(value - nearest) > integerTolerance)
1002                unsatisfied++;
1003        }
1004    }
1005    int way = branchingObject->way();
1006    way = - way; // because after branch so moved on
1007    double value = branchingObject->value();
1008    CbcObjectUpdateData newData (this, way,
1009                                 change, iStatus,
1010                                 originalUnsatisfied - unsatisfied, value);
1011    newData.originalObjective_ = originalValue;
1012    // Solvers know about direction
1013    double direction = solver->getObjSense();
1014    solver->getDblParam(OsiDualObjectiveLimit, newData.cutoff_);
1015    newData.cutoff_ *= direction;
1016    return newData;
1017}
1018// Just update using feasible branches and keep count of infeasible
1019#undef INFEAS
1020// Update object by CbcObjectUpdateData
1021void
1022CbcSimpleIntegerDynamicPseudoCost::updateInformation(const CbcObjectUpdateData & data)
1023{
1024    bool feasible = data.status_ != 1;
1025    int way = data.way_;
1026    double value = data.branchingValue_;
1027    double change = data.change_;
1028#ifdef COIN_DEVELOP
1029    History hist;
1030    hist.where_ = 'U'; // need to tell if hot
1031#endif
1032    double movement = 0.0;
1033    if (way < 0) {
1034        // down
1035        movement = value - floor(value);
1036        if (feasible) {
1037#ifdef COIN_DEVELOP
1038            hist.status_ = 'D';
1039#endif
1040            movement = CoinMax(movement, MINIMUM_MOVEMENT);
1041            //printf("(down change %g value down %g ",change,movement);
1042            incrementNumberTimesDown();
1043            addToSumDownChange(1.0e-30 + movement);
1044            addToSumDownDecrease(data.intDecrease_);
1045#if TYPE2==0
1046            addToSumDownCost(change / (1.0e-30 + movement));
1047            setDownDynamicPseudoCost(sumDownCost() / static_cast<double>( numberTimesDown()));
1048#elif TYPE2==1
1049            addToSumDownCost(change);
1050            setDownDynamicPseudoCost(sumDownCost() / sumDownChange());
1051#elif TYPE2==2
1052            addToSumDownCost(change*TYPERATIO + (1.0 - TYPERATIO)*change / (1.0e-30 + movement));
1053            setDownDynamicPseudoCost(sumDownCost()*(TYPERATIO / sumDownChange() + (1.0 - TYPERATIO) / (double) numberTimesDown()));
1054#endif
1055        } else {
1056#ifdef COIN_DEVELOP
1057            hist.status_ = 'd';
1058#endif
1059            //printf("(down infeasible value down %g ",change,movement);
1060            incrementNumberTimesDown();
1061            incrementNumberTimesDownInfeasible();
1062#if INFEAS==2
1063            double distanceToCutoff = 0.0;
1064            double objectiveValue = model->getCurrentMinimizationObjValue();
1065            distanceToCutoff =  model->getCutoff()  - originalValue;
1066            if (distanceToCutoff < 1.0e20)
1067                change = distanceToCutoff * 2.0;
1068            else
1069                change = downDynamicPseudoCost() * movement * 10.0;
1070            change = CoinMax(1.0e-12 * (1.0 + fabs(originalValue)), change);
1071            addToSumDownChange(1.0e-30 + movement);
1072            addToSumDownDecrease(data.intDecrease_);
1073#if TYPE2==0
1074            addToSumDownCost(change / (1.0e-30 + movement));
1075            setDownDynamicPseudoCost(sumDownCost() / (double) numberTimesDown());
1076#elif TYPE2==1
1077            addToSumDownCost(change);
1078            setDownDynamicPseudoCost(sumDownCost() / sumDownChange());
1079#elif TYPE2==2
1080            addToSumDownCost(change*TYPERATIO + (1.0 - TYPERATIO)*change / (1.0e-30 + movement));
1081            setDownDynamicPseudoCost(sumDownCost()*(TYPERATIO / sumDownChange() + (1.0 - TYPERATIO) / (double) numberTimesDown()));
1082#endif
1083#endif
1084        }
1085#if INFEAS==1
1086        double sum = sumDownCost_;
1087        int number = numberTimesDown_;
1088        double originalValue = data.originalObjective_;
1089        assert (originalValue != COIN_DBL_MAX);
1090        double distanceToCutoff =  data.cutoff_  - originalValue;
1091        if (distanceToCutoff > 1.0e20)
1092            distanceToCutoff = 10.0 + fabs(originalValue);
1093        sum += numberTimesDownInfeasible_ * CoinMax(distanceToCutoff, 1.0e-12 * (1.0 + fabs(originalValue)));
1094        setDownDynamicPseudoCost(sum / static_cast<double> (number));
1095#endif
1096    } else {
1097        // up
1098        movement = ceil(value) - value;
1099        if (feasible) {
1100#ifdef COIN_DEVELOP
1101            hist.status_ = 'U';
1102#endif
1103            movement = CoinMax(movement, MINIMUM_MOVEMENT);
1104            //printf("(up change %g value down %g ",change,movement);
1105            incrementNumberTimesUp();
1106            addToSumUpChange(1.0e-30 + movement);
1107            addToSumUpDecrease(data.intDecrease_);
1108#if TYPE2==0
1109            addToSumUpCost(change / (1.0e-30 + movement));
1110            setUpDynamicPseudoCost(sumUpCost() / static_cast<double> (numberTimesUp()));
1111#elif TYPE2==1
1112            addToSumUpCost(change);
1113            setUpDynamicPseudoCost(sumUpCost() / sumUpChange());
1114#elif TYPE2==2
1115            addToSumUpCost(change*TYPERATIO + (1.0 - TYPERATIO)*change / (1.0e-30 + movement));
1116            setUpDynamicPseudoCost(sumUpCost()*(TYPERATIO / sumUpChange() + (1.0 - TYPERATIO) / (double) numberTimesUp()));
1117#endif
1118        } else {
1119#ifdef COIN_DEVELOP
1120            hist.status_ = 'u';
1121#endif
1122            //printf("(up infeasible value down %g ",change,movement);
1123            incrementNumberTimesUp();
1124            incrementNumberTimesUpInfeasible();
1125#if INFEAS==2
1126            double distanceToCutoff = 0.0;
1127            double objectiveValue = model->getCurrentMinimizationObjValue();
1128            distanceToCutoff =  model->getCutoff()  - originalValue;
1129            if (distanceToCutoff < 1.0e20)
1130                change = distanceToCutoff * 2.0;
1131            else
1132                change = upDynamicPseudoCost() * movement * 10.0;
1133            change = CoinMax(1.0e-12 * (1.0 + fabs(originalValue)), change);
1134            addToSumUpChange(1.0e-30 + movement);
1135            addToSumUpDecrease(data.intDecrease_);
1136#if TYPE2==0
1137            addToSumUpCost(change / (1.0e-30 + movement));
1138            setUpDynamicPseudoCost(sumUpCost() / (double) numberTimesUp());
1139#elif TYPE2==1
1140            addToSumUpCost(change);
1141            setUpDynamicPseudoCost(sumUpCost() / sumUpChange());
1142#elif TYPE2==2
1143            addToSumUpCost(change*TYPERATIO + (1.0 - TYPERATIO)*change / (1.0e-30 + movement));
1144            setUpDynamicPseudoCost(sumUpCost()*(TYPERATIO / sumUpChange() + (1.0 - TYPERATIO) / (double) numberTimesUp()));
1145#endif
1146#endif
1147        }
1148#if INFEAS==1
1149        double sum = sumUpCost_;
1150        int number = numberTimesUp_;
1151        double originalValue = data.originalObjective_;
1152        assert (originalValue != COIN_DBL_MAX);
1153        double distanceToCutoff =  data.cutoff_  - originalValue;
1154        if (distanceToCutoff > 1.0e20)
1155            distanceToCutoff = 10.0 + fabs(originalValue);
1156        sum += numberTimesUpInfeasible_ * CoinMax(distanceToCutoff, 1.0e-12 * (1.0 + fabs(originalValue)));
1157        setUpDynamicPseudoCost(sum / static_cast<double> (number));
1158#endif
1159    }
1160    if (data.way_ < 0)
1161        assert (numberTimesDown_ > 0);
1162    else
1163        assert (numberTimesUp_ > 0);
1164    assert (downDynamicPseudoCost_ >= 0.0 && downDynamicPseudoCost_ < 1.0e100);
1165    downDynamicPseudoCost_ = CoinMax(1.0e-10, downDynamicPseudoCost_);
1166    assert (upDynamicPseudoCost_ >= 0.0 && upDynamicPseudoCost_ < 1.0e100);
1167    upDynamicPseudoCost_ = CoinMax(1.0e-10, upDynamicPseudoCost_);
1168#ifdef COIN_DEVELOP
1169    hist.sequence_ = columnNumber_;
1170    hist.numberUp_ = numberTimesUp_;
1171    hist.numberUpInf_ = numberTimesUpInfeasible_;
1172    hist.sumUp_ = sumUpCost_;
1173    hist.upEst_ = change;
1174    hist.numberDown_ = numberTimesDown_;
1175    hist.numberDownInf_ = numberTimesDownInfeasible_;
1176    hist.sumDown_ = sumDownCost_;
1177    hist.downEst_ = movement;
1178    addRecord(hist);
1179#endif
1180    //print(1,0.5);
1181    assert (downDynamicPseudoCost_ > 1.0e-40 && upDynamicPseudoCost_ > 1.0e-40);
1182#if MOD_SHADOW>1
1183    if (upShadowPrice_ > 0.0 && numberTimesDown_ >= numberBeforeTrust_
1184            && numberTimesUp_ >= numberBeforeTrust_) {
1185        // Set negative
1186        upShadowPrice_ = -upShadowPrice_;
1187        assert (downShadowPrice_ > 0.0);
1188        downShadowPrice_ = - downShadowPrice_;
1189    }
1190#endif
1191}
1192// Updates stuff like pseudocosts after mini branch and bound
1193void
1194CbcSimpleIntegerDynamicPseudoCost::updateAfterMini(int numberDown, int numberDownInfeasible,
1195        double sumDown, int numberUp,
1196        int numberUpInfeasible, double sumUp)
1197{
1198    numberTimesDown_ = numberDown;
1199    numberTimesDownInfeasible_ = numberDownInfeasible;
1200    sumDownCost_ = sumDown;
1201    numberTimesUp_ = numberUp;
1202    numberTimesUpInfeasible_ = numberUpInfeasible;
1203    sumUpCost_ = sumUp;
1204    if (numberTimesDown_ > 0) {
1205        setDownDynamicPseudoCost(sumDownCost_ / static_cast<double> (numberTimesDown_));
1206        assert (downDynamicPseudoCost_ > 0.0 && downDynamicPseudoCost_ < 1.0e50);
1207    }
1208    if (numberTimesUp_ > 0) {
1209        setUpDynamicPseudoCost(sumUpCost_ / static_cast<double> (numberTimesUp_));
1210        assert (upDynamicPseudoCost_ > 0.0 && upDynamicPseudoCost_ < 1.0e50);
1211    }
1212    assert (downDynamicPseudoCost_ > 1.0e-40 && upDynamicPseudoCost_ > 1.0e-40);
1213}
1214// Pass in probing information
1215void
1216CbcSimpleIntegerDynamicPseudoCost::setProbingInformation(int fixedDown, int fixedUp)
1217{
1218    numberTimesProbingTotal_++;
1219    numberTimesDownLocalFixed_ = fixedDown;
1220    numberTimesDownTotalFixed_ += fixedDown;
1221    numberTimesUpLocalFixed_ = fixedUp;
1222    numberTimesUpTotalFixed_ += fixedUp;
1223}
1224// Print
1225void
1226CbcSimpleIntegerDynamicPseudoCost::print(int type, double value) const
1227{
1228    if (!type) {
1229        double meanDown = 0.0;
1230        double devDown = 0.0;
1231        if (numberTimesDown_) {
1232            meanDown = sumDownCost_ / static_cast<double> (numberTimesDown_);
1233            devDown = meanDown * meanDown - 2.0 * meanDown * sumDownCost_;
1234            if (devDown >= 0.0)
1235                devDown = sqrt(devDown);
1236        }
1237        double meanUp = 0.0;
1238        double devUp = 0.0;
1239        if (numberTimesUp_) {
1240            meanUp = sumUpCost_ / static_cast<double> (numberTimesUp_);
1241            devUp = meanUp * meanUp - 2.0 * meanUp * sumUpCost_;
1242            if (devUp >= 0.0)
1243                devUp = sqrt(devUp);
1244        }
1245        printf("%d down %d times (%d inf) mean %g (dev %g) up %d times (%d inf) mean %g (dev %g)\n",
1246               columnNumber_,
1247               numberTimesDown_, numberTimesDownInfeasible_, meanDown, devDown,
1248               numberTimesUp_, numberTimesUpInfeasible_, meanUp, devUp);
1249    } else {
1250        const double * upper = model_->getCbcColUpper();
1251        double integerTolerance =
1252            model_->getDblParam(CbcModel::CbcIntegerTolerance);
1253        double below = floor(value + integerTolerance);
1254        double above = below + 1.0;
1255        if (above > upper[columnNumber_]) {
1256            above = below;
1257            below = above - 1;
1258        }
1259        double objectiveValue = model_->getCurrentMinimizationObjValue();
1260        double distanceToCutoff =  model_->getCutoff() - objectiveValue;
1261        if (distanceToCutoff < 1.0e20)
1262            distanceToCutoff *= 10.0;
1263        else
1264            distanceToCutoff = 1.0e2 + fabs(objectiveValue);
1265        distanceToCutoff = CoinMax(distanceToCutoff, 1.0e-12 * (1.0 + fabs(objectiveValue)));
1266        double sum;
1267        int number;
1268        double downCost = CoinMax(value - below, 0.0);
1269        double downCost0 = downCost * downDynamicPseudoCost_;
1270        sum = sumDownCost();
1271        number = numberTimesDown();
1272        sum += numberTimesDownInfeasible() * (distanceToCutoff / (downCost + 1.0e-12));
1273        if (number > 0)
1274            downCost *= sum / static_cast<double> (number);
1275        else
1276            downCost  *=  downDynamicPseudoCost_;
1277        double upCost = CoinMax((above - value), 0.0);
1278        double upCost0 = upCost * upDynamicPseudoCost_;
1279        sum = sumUpCost();
1280        number = numberTimesUp();
1281        sum += numberTimesUpInfeasible() * (distanceToCutoff / (upCost + 1.0e-12));
1282        if (number > 0)
1283            upCost *= sum / static_cast<double> (number);
1284        else
1285            upCost  *=  upDynamicPseudoCost_;
1286        printf("%d down %d times %g (est %g)  up %d times %g (est %g)\n",
1287               columnNumber_,
1288               numberTimesDown_, downCost, downCost0,
1289               numberTimesUp_, upCost, upCost0);
1290    }
1291}
1292
1293//##############################################################################
1294
1295// Default Constructor
1296CbcIntegerPseudoCostBranchingObject::CbcIntegerPseudoCostBranchingObject()
1297        : CbcIntegerBranchingObject()
1298{
1299    changeInGuessed_ = 1.0e-5;
1300}
1301
1302// Useful constructor
1303CbcIntegerPseudoCostBranchingObject::CbcIntegerPseudoCostBranchingObject (CbcModel * model,
1304        int variable, int way , double value)
1305        : CbcIntegerBranchingObject(model, variable, way, value)
1306{
1307}
1308// Useful constructor for fixing
1309CbcIntegerPseudoCostBranchingObject::CbcIntegerPseudoCostBranchingObject (CbcModel * model,
1310        int variable, int way,
1311        double lowerValue,
1312        double /*upperValue*/)
1313        : CbcIntegerBranchingObject(model, variable, way, lowerValue)
1314{
1315    changeInGuessed_ = 1.0e100;
1316}
1317
1318
1319// Copy constructor
1320CbcIntegerPseudoCostBranchingObject::CbcIntegerPseudoCostBranchingObject (
1321    const CbcIntegerPseudoCostBranchingObject & rhs)
1322        : CbcIntegerBranchingObject(rhs)
1323{
1324    changeInGuessed_ = rhs.changeInGuessed_;
1325}
1326
1327// Assignment operator
1328CbcIntegerPseudoCostBranchingObject &
1329CbcIntegerPseudoCostBranchingObject::operator=( const CbcIntegerPseudoCostBranchingObject & rhs)
1330{
1331    if (this != &rhs) {
1332        CbcIntegerBranchingObject::operator=(rhs);
1333        changeInGuessed_ = rhs.changeInGuessed_;
1334    }
1335    return *this;
1336}
1337CbcBranchingObject *
1338CbcIntegerPseudoCostBranchingObject::clone() const
1339{
1340    return (new CbcIntegerPseudoCostBranchingObject(*this));
1341}
1342
1343
1344// Destructor
1345CbcIntegerPseudoCostBranchingObject::~CbcIntegerPseudoCostBranchingObject ()
1346{
1347}
1348
1349/*
1350  Perform a branch by adjusting the bounds of the specified variable. Note
1351  that each arm of the branch advances the object to the next arm by
1352  advancing the value of way_.
1353
1354  Providing new values for the variable's lower and upper bounds for each
1355  branching direction gives a little bit of additional flexibility and will
1356  be easily extensible to multi-way branching.
1357  Returns change in guessed objective on next branch
1358*/
1359double
1360CbcIntegerPseudoCostBranchingObject::branch()
1361{
1362    CbcIntegerBranchingObject::branch();
1363    return changeInGuessed_;
1364}
1365
1366/** Compare the \c this with \c brObj. \c this and \c brObj must be os the
1367    same type and must have the same original object, but they may have
1368    different feasible regions.
1369    Return the appropriate CbcRangeCompare value (first argument being the
1370    sub/superset if that's the case). In case of overlap (and if \c
1371    replaceIfOverlap is true) replace the current branching object with one
1372    whose feasible region is the overlap.
1373*/
1374CbcRangeCompare
1375CbcIntegerPseudoCostBranchingObject::compareBranchingObject
1376(const CbcBranchingObject* brObj, const bool replaceIfOverlap)
1377{
1378    const CbcIntegerPseudoCostBranchingObject* br =
1379        dynamic_cast<const CbcIntegerPseudoCostBranchingObject*>(brObj);
1380    assert(br);
1381    double* thisBd = way_ < 0 ? down_ : up_;
1382    const double* otherBd = br->way_ < 0 ? br->down_ : br->up_;
1383    return CbcCompareRanges(thisBd, otherBd, replaceIfOverlap);
1384}
Note: See TracBrowser for help on using the repository browser.