source: branches/sandbox/Cbc/src/CbcSimpleIntegerDynamicPseudoCost.cpp @ 1308

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

Broke up CbcBranchDynamic? and CbcBranchLotsize?.cpp.
Updated spreadsheets.

File size: 49.0 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    double value = solution[columnNumber_];
574    value = CoinMax(value, lower[columnNumber_]);
575    value = CoinMin(value, upper[columnNumber_]);
576    /*printf("%d %g %g %g %g\n",columnNumber_,value,lower[columnNumber_],
577      solution[columnNumber_],upper[columnNumber_]);*/
578    double nearest = floor(value + 0.5);
579    double integerTolerance =
580        model_->getDblParam(CbcModel::CbcIntegerTolerance);
581    double below = floor(value + integerTolerance);
582    double above = below + 1.0;
583    if (above > upper[columnNumber_]) {
584        above = below;
585        below = above - 1;
586    }
587#if INFEAS==1
588    double distanceToCutoff = 0.0;
589    double objectiveValue = model_->getCurrentMinimizationObjValue();
590    distanceToCutoff =  model_->getCutoff()  - objectiveValue;
591    if (distanceToCutoff < 1.0e20)
592        distanceToCutoff *= 10.0;
593    else
594        distanceToCutoff = 1.0e2 + fabs(objectiveValue);
595    distanceToCutoff = CoinMax(distanceToCutoff, 1.0e-12 * (1.0 + fabs(objectiveValue)));
596#endif
597    double sum;
598    double number;
599    double downCost = CoinMax(value - below, 0.0);
600#if TYPE2==0
601    sum = sumDownCost_;
602    number = numberTimesDown_;
603#if INFEAS==1
604    sum += numberTimesDownInfeasible_ * CoinMax(distanceToCutoff / (downCost + 1.0e-12), sumDownCost_);
605#endif
606#elif TYPE2==1
607    sum = sumDownCost_;
608    number = sumDownChange_;
609#if INFEAS==1
610    sum += numberTimesDownInfeasible_ * CoinMax(distanceToCutoff / (downCost + 1.0e-12), sumDownCost_);
611#endif
612#elif TYPE2==2
613    abort();
614#if INFEAS==1
615    sum += numberTimesDownInfeasible_ * (distanceToCutoff / (downCost + 1.0e-12));
616#endif
617#endif
618#if MOD_SHADOW>0
619    if (!downShadowPrice_) {
620        if (number > 0.0)
621            downCost *= sum / number;
622        else
623            downCost  *=  downDynamicPseudoCost_;
624    } else if (downShadowPrice_ > 0.0) {
625        downCost *= downShadowPrice_;
626    } else {
627        downCost *= (downDynamicPseudoCost_ - downShadowPrice_);
628    }
629#else
630    if (downShadowPrice_ <= 0.0) {
631        if (number > 0.0)
632            downCost *= sum / number;
633        else
634            downCost  *=  downDynamicPseudoCost_;
635    } else {
636        downCost *= downShadowPrice_;
637    }
638#endif
639    double upCost = CoinMax((above - value), 0.0);
640#if TYPE2==0
641    sum = sumUpCost_;
642    number = numberTimesUp_;
643#if INFEAS==1
644    sum += numberTimesUpInfeasible_ * CoinMax(distanceToCutoff / (upCost + 1.0e-12), sumUpCost_);
645#endif
646#elif TYPE2==1
647    sum = sumUpCost_;
648    number = sumUpChange_;
649#if INFEAS==1
650    sum += numberTimesUpInfeasible_ * CoinMax(distanceToCutoff / (upCost + 1.0e-12), sumUpCost_);
651#endif
652#elif TYPE2==1
653    abort();
654#if INFEAS==1
655    sum += numberTimesUpInfeasible_ * (distanceToCutoff / (upCost + 1.0e-12));
656#endif
657#endif
658#if MOD_SHADOW>0
659    if (!upShadowPrice_) {
660        if (number > 0.0)
661            upCost *= sum / number;
662        else
663            upCost  *=  upDynamicPseudoCost_;
664    } else if (upShadowPrice_ > 0.0) {
665        upCost *= upShadowPrice_;
666    } else {
667        upCost *= (upDynamicPseudoCost_ - upShadowPrice_);
668    }
669#else
670    if (upShadowPrice_ <= 0.0) {
671        if (number > 0.0)
672            upCost *= sum / number;
673        else
674            upCost  *=  upDynamicPseudoCost_;
675    } else {
676        upCost *= upShadowPrice_;
677    }
678#endif
679    if (downCost >= upCost)
680        preferredWay = 1;
681    else
682        preferredWay = -1;
683    // See if up down choice set
684    if (upDownSeparator_ > 0.0) {
685        preferredWay = (value - below >= upDownSeparator_) ? 1 : -1;
686    }
687#ifdef FUNNY_BRANCHING
688    if (fabs(value - nearest) > integerTolerance) {
689        double ratio = (100.0 + lastUpDecrease_) / (100.0 + lastDownDecrease_);
690        downCost *= ratio;
691        upCost *= ratio;
692        if ((lastUpDecrease_ % 100) == -1)
693            printf("col %d total %d djtimes %d\n",
694                   columnNumber_, lastDownDecrease_, lastUpDecrease_);
695    }
696#endif
697    if (preferredWay_)
698        preferredWay = preferredWay_;
699    if (info->hotstartSolution_) {
700        double targetValue = info->hotstartSolution_[columnNumber_];
701        if (value > targetValue)
702            preferredWay = -1;
703        else
704            preferredWay = 1;
705    }
706    if (fabs(value - nearest) <= integerTolerance) {
707        if (priority_ != -999)
708            return 0.0;
709        else
710            return 1.0e-13;
711    } else {
712        int stateOfSearch = model_->stateOfSearch() % 10;
713        double returnValue = 0.0;
714        double minValue = CoinMin(downCost, upCost);
715        double maxValue = CoinMax(downCost, upCost);
716        char where;
717        // was <= 10
718        //if (stateOfSearch<=1||model_->currentNode()->depth()<=-10 /* was ||maxValue>0.2*distanceToCutoff*/) {
719        if (stateOfSearch <= 2) {
720            // no branching solution
721            where = 'i';
722            returnValue = WEIGHT_BEFORE * minValue + (1.0 - WEIGHT_BEFORE) * maxValue;
723            if (0) {
724                double sum;
725                int number;
726                double downCost2 = CoinMax(value - below, 0.0);
727                sum = sumDownCost_;
728                number = numberTimesDown_;
729                if (number > 0)
730                    downCost2 *= sum / static_cast<double> (number);
731                else
732                    downCost2  *=  downDynamicPseudoCost_;
733                double upCost2 = CoinMax((above - value), 0.0);
734                sum = sumUpCost_;
735                number = numberTimesUp_;
736                if (number > 0)
737                    upCost2 *= sum / static_cast<double> (number);
738                else
739                    upCost2  *=  upDynamicPseudoCost_;
740                double minValue2 = CoinMin(downCost2, upCost2);
741                double maxValue2 = CoinMax(downCost2, upCost2);
742                printf("%d value %g downC %g upC %g minV %g maxV %g downC2 %g upC2 %g minV2 %g maxV2 %g\n",
743                       columnNumber_, value, downCost, upCost, minValue, maxValue,
744                       downCost2, upCost2, minValue2, maxValue2);
745            }
746        } else {
747            // some solution
748            where = 'I';
749#ifndef WEIGHT_PRODUCT
750            returnValue = WEIGHT_AFTER * minValue + (1.0 - WEIGHT_AFTER) * maxValue;
751#else
752            double minProductWeight = model_->getDblParam(CbcModel::CbcSmallChange);
753            returnValue = CoinMax(minValue, minProductWeight) * CoinMax(maxValue, minProductWeight);
754            //returnValue += minProductWeight*minValue;
755#endif
756        }
757        if (numberTimesUp_ < numberBeforeTrust_ ||
758                numberTimesDown_ < numberBeforeTrust_) {
759            //if (returnValue<1.0e10)
760            //returnValue += 1.0e12;
761            //else
762            returnValue *= 1.0e3;
763            if (!numberTimesUp_ && !numberTimesDown_)
764                returnValue *= 1.0e10;
765        }
766        //if (fabs(value-0.5)<1.0e-5) {
767        //returnValue = 3.0*returnValue + 0.2;
768        //} else if (value>0.9) {
769        //returnValue = 2.0*returnValue + 0.1;
770        //}
771        if (method_ == 1) {
772            // probing
773            // average
774            double up = 1.0e-15;
775            double down = 1.0e-15;
776            if (numberTimesProbingTotal_) {
777                up += numberTimesUpTotalFixed_ / static_cast<double> (numberTimesProbingTotal_);
778                down += numberTimesDownTotalFixed_ / static_cast<double> (numberTimesProbingTotal_);
779            }
780            returnValue = 1 + 10.0 * CoinMin(numberTimesDownLocalFixed_, numberTimesUpLocalFixed_) +
781                          CoinMin(down, up);
782            returnValue *= 1.0e-3;
783        }
784#ifdef COIN_DEVELOP
785        History hist;
786        hist.where_ = where;
787        hist.status_ = ' ';
788        hist.sequence_ = columnNumber_;
789        hist.numberUp_ = numberTimesUp_;
790        hist.numberUpInf_ = numberTimesUpInfeasible_;
791        hist.sumUp_ = sumUpCost_;
792        hist.upEst_ = upCost;
793        hist.numberDown_ = numberTimesDown_;
794        hist.numberDownInf_ = numberTimesDownInfeasible_;
795        hist.sumDown_ = sumDownCost_;
796        hist.downEst_ = downCost;
797        if (stateOfSearch)
798            addRecord(hist);
799#endif
800        return CoinMax(returnValue, 1.0e-15);
801    }
802}
803// Creates a branching object
804CbcBranchingObject *
805CbcSimpleIntegerDynamicPseudoCost::createCbcBranch(OsiSolverInterface * /*solver*/,
806        const OsiBranchingInformation * info, int way)
807{
808    double value = info->solution_[columnNumber_];
809    value = CoinMax(value, info->lower_[columnNumber_]);
810    value = CoinMin(value, info->upper_[columnNumber_]);
811    assert (info->upper_[columnNumber_] > info->lower_[columnNumber_]);
812    if (!info->hotstartSolution_ && priority_ != -999) {
813#ifndef NDEBUG
814        double nearest = floor(value + 0.5);
815        assert (fabs(value - nearest) > info->integerTolerance_);
816#endif
817    } else if (info->hotstartSolution_) {
818        double targetValue = info->hotstartSolution_[columnNumber_];
819        if (way > 0)
820            value = targetValue - 0.1;
821        else
822            value = targetValue + 0.1;
823    } else {
824        if (value <= info->lower_[columnNumber_])
825            value += 0.1;
826        else if (value >= info->upper_[columnNumber_])
827            value -= 0.1;
828    }
829    assert (value >= info->lower_[columnNumber_] &&
830            value <= info->upper_[columnNumber_]);
831    CbcDynamicPseudoCostBranchingObject * newObject =
832        new CbcDynamicPseudoCostBranchingObject(model_, columnNumber_, way,
833                                                value, this);
834    double up =  upDynamicPseudoCost_ * (ceil(value) - value);
835    double down =  downDynamicPseudoCost_ * (value - floor(value));
836    double changeInGuessed = up - down;
837    if (way > 0)
838        changeInGuessed = - changeInGuessed;
839    changeInGuessed = CoinMax(0.0, changeInGuessed);
840    //if (way>0)
841    //changeInGuessed += 1.0e8; // bias to stay up
842    newObject->setChangeInGuessed(changeInGuessed);
843    newObject->setOriginalObject(this);
844    return newObject;
845}
846
847// Return "up" estimate
848double
849CbcSimpleIntegerDynamicPseudoCost::upEstimate() const
850{
851    const double * solution = model_->testSolution();
852    const double * lower = model_->getCbcColLower();
853    const double * upper = model_->getCbcColUpper();
854    double value = solution[columnNumber_];
855    value = CoinMax(value, lower[columnNumber_]);
856    value = CoinMin(value, upper[columnNumber_]);
857    if (upper[columnNumber_] == lower[columnNumber_]) {
858        // fixed
859        return 0.0;
860    }
861    double integerTolerance =
862        model_->getDblParam(CbcModel::CbcIntegerTolerance);
863    double below = floor(value + integerTolerance);
864    double above = below + 1.0;
865    if (above > upper[columnNumber_]) {
866        above = below;
867        below = above - 1;
868    }
869    double upCost = CoinMax((above - value) * upDynamicPseudoCost_, 0.0);
870    return upCost;
871}
872// Return "down" estimate
873double
874CbcSimpleIntegerDynamicPseudoCost::downEstimate() const
875{
876    const double * solution = model_->testSolution();
877    const double * lower = model_->getCbcColLower();
878    const double * upper = model_->getCbcColUpper();
879    double value = solution[columnNumber_];
880    value = CoinMax(value, lower[columnNumber_]);
881    value = CoinMin(value, upper[columnNumber_]);
882    if (upper[columnNumber_] == lower[columnNumber_]) {
883        // fixed
884        return 0.0;
885    }
886    double integerTolerance =
887        model_->getDblParam(CbcModel::CbcIntegerTolerance);
888    double below = floor(value + integerTolerance);
889    double above = below + 1.0;
890    if (above > upper[columnNumber_]) {
891        above = below;
892        below = above - 1;
893    }
894    double downCost = CoinMax((value - below) * downDynamicPseudoCost_, 0.0);
895    return downCost;
896}
897// Set down pseudo cost
898void
899CbcSimpleIntegerDynamicPseudoCost::setDownDynamicPseudoCost(double value)
900{
901#ifdef TRACE_ONE
902    double oldDown = sumDownCost_;
903#endif
904    downDynamicPseudoCost_ = value;
905    sumDownCost_ = CoinMax(sumDownCost_, value * numberTimesDown_);
906#ifdef TRACE_ONE
907    if (columnNumber_ == TRACE_ONE) {
908        double down = downDynamicPseudoCost_ * numberTimesDown_;
909        printf("For %d sumDown %g (%d), inf (%d) - pseudo %g - sumDown was %g -> %g\n",
910               TRACE_ONE, down, numberTimesDown_,
911               numberTimesDownInfeasible_, downDynamicPseudoCost_,
912               oldDown, sumDownCost_);
913    }
914#endif
915}
916// Modify down pseudo cost in a slightly different way
917void
918CbcSimpleIntegerDynamicPseudoCost::updateDownDynamicPseudoCost(double value)
919{
920    sumDownCost_ += value;
921    numberTimesDown_++;
922    downDynamicPseudoCost_ = sumDownCost_ / static_cast<double>(numberTimesDown_);
923}
924// Set up pseudo cost
925void
926CbcSimpleIntegerDynamicPseudoCost::setUpDynamicPseudoCost(double value)
927{
928#ifdef TRACE_ONE
929    double oldUp = sumUpCost_;
930#endif
931    upDynamicPseudoCost_ = value;
932    sumUpCost_ = CoinMax(sumUpCost_, value * numberTimesUp_);
933#ifdef TRACE_ONE
934    if (columnNumber_ == TRACE_ONE) {
935        double up = upDynamicPseudoCost_ * numberTimesUp_;
936        printf("For %d sumUp %g (%d), inf (%d) - pseudo %g - sumUp was %g -> %g\n",
937               TRACE_ONE, up, numberTimesUp_,
938               numberTimesUpInfeasible_, upDynamicPseudoCost_,
939               oldUp, sumUpCost_);
940    }
941#endif
942}
943// Modify up pseudo cost in a slightly different way
944void
945CbcSimpleIntegerDynamicPseudoCost::updateUpDynamicPseudoCost(double value)
946{
947    sumUpCost_ += value;
948    numberTimesUp_++;
949    upDynamicPseudoCost_ = sumUpCost_ / static_cast<double>(numberTimesUp_);
950}
951/* Pass in information on branch just done and create CbcObjectUpdateData instance.
952   If object does not need data then backward pointer will be NULL.
953   Assumes can get information from solver */
954CbcObjectUpdateData
955CbcSimpleIntegerDynamicPseudoCost::createUpdateInformation(const OsiSolverInterface * solver,
956        const CbcNode * node,
957        const CbcBranchingObject * branchingObject)
958{
959    double originalValue = node->objectiveValue();
960    int originalUnsatisfied = node->numberUnsatisfied();
961    double objectiveValue = solver->getObjValue() * solver->getObjSense();
962    int unsatisfied = 0;
963    int i;
964    //might be base model - doesn't matter
965    int numberIntegers = model_->numberIntegers();;
966    const double * solution = solver->getColSolution();
967    //const double * lower = solver->getColLower();
968    //const double * upper = solver->getColUpper();
969    double change = CoinMax(0.0, objectiveValue - originalValue);
970    int iStatus;
971    if (solver->isProvenOptimal())
972        iStatus = 0; // optimal
973    else if (solver->isIterationLimitReached()
974             && !solver->isDualObjectiveLimitReached())
975        iStatus = 2; // unknown
976    else
977        iStatus = 1; // infeasible
978
979    bool feasible = iStatus != 1;
980    if (feasible) {
981        double integerTolerance =
982            model_->getDblParam(CbcModel::CbcIntegerTolerance);
983        const int * integerVariable = model_->integerVariable();
984        for (i = 0; i < numberIntegers; i++) {
985            int j = integerVariable[i];
986            double value = solution[j];
987            double nearest = floor(value + 0.5);
988            if (fabs(value - nearest) > integerTolerance)
989                unsatisfied++;
990        }
991    }
992    int way = branchingObject->way();
993    way = - way; // because after branch so moved on
994    double value = branchingObject->value();
995    CbcObjectUpdateData newData (this, way,
996                                 change, iStatus,
997                                 originalUnsatisfied - unsatisfied, value);
998    newData.originalObjective_ = originalValue;
999    // Solvers know about direction
1000    double direction = solver->getObjSense();
1001    solver->getDblParam(OsiDualObjectiveLimit, newData.cutoff_);
1002    newData.cutoff_ *= direction;
1003    return newData;
1004}
1005// Just update using feasible branches and keep count of infeasible
1006#undef INFEAS
1007// Update object by CbcObjectUpdateData
1008void
1009CbcSimpleIntegerDynamicPseudoCost::updateInformation(const CbcObjectUpdateData & data)
1010{
1011    bool feasible = data.status_ != 1;
1012    int way = data.way_;
1013    double value = data.branchingValue_;
1014    double change = data.change_;
1015#ifdef COIN_DEVELOP
1016    History hist;
1017    hist.where_ = 'U'; // need to tell if hot
1018#endif
1019    double movement = 0.0;
1020    if (way < 0) {
1021        // down
1022        movement = value - floor(value);
1023        if (feasible) {
1024#ifdef COIN_DEVELOP
1025            hist.status_ = 'D';
1026#endif
1027            movement = CoinMax(movement, MINIMUM_MOVEMENT);
1028            //printf("(down change %g value down %g ",change,movement);
1029            incrementNumberTimesDown();
1030            addToSumDownChange(1.0e-30 + movement);
1031            addToSumDownDecrease(data.intDecrease_);
1032#if TYPE2==0
1033            addToSumDownCost(change / (1.0e-30 + movement));
1034            setDownDynamicPseudoCost(sumDownCost() / static_cast<double>( numberTimesDown()));
1035#elif TYPE2==1
1036            addToSumDownCost(change);
1037            setDownDynamicPseudoCost(sumDownCost() / sumDownChange());
1038#elif TYPE2==2
1039            addToSumDownCost(change*TYPERATIO + (1.0 - TYPERATIO)*change / (1.0e-30 + movement));
1040            setDownDynamicPseudoCost(sumDownCost()*(TYPERATIO / sumDownChange() + (1.0 - TYPERATIO) / (double) numberTimesDown()));
1041#endif
1042        } else {
1043#ifdef COIN_DEVELOP
1044            hist.status_ = 'd';
1045#endif
1046            //printf("(down infeasible value down %g ",change,movement);
1047            incrementNumberTimesDown();
1048            incrementNumberTimesDownInfeasible();
1049#if INFEAS==2
1050            double distanceToCutoff = 0.0;
1051            double objectiveValue = model->getCurrentMinimizationObjValue();
1052            distanceToCutoff =  model->getCutoff()  - originalValue;
1053            if (distanceToCutoff < 1.0e20)
1054                change = distanceToCutoff * 2.0;
1055            else
1056                change = downDynamicPseudoCost() * movement * 10.0;
1057            change = CoinMax(1.0e-12 * (1.0 + fabs(originalValue)), change);
1058            addToSumDownChange(1.0e-30 + movement);
1059            addToSumDownDecrease(data.intDecrease_);
1060#if TYPE2==0
1061            addToSumDownCost(change / (1.0e-30 + movement));
1062            setDownDynamicPseudoCost(sumDownCost() / (double) numberTimesDown());
1063#elif TYPE2==1
1064            addToSumDownCost(change);
1065            setDownDynamicPseudoCost(sumDownCost() / sumDownChange());
1066#elif TYPE2==2
1067            addToSumDownCost(change*TYPERATIO + (1.0 - TYPERATIO)*change / (1.0e-30 + movement));
1068            setDownDynamicPseudoCost(sumDownCost()*(TYPERATIO / sumDownChange() + (1.0 - TYPERATIO) / (double) numberTimesDown()));
1069#endif
1070#endif
1071        }
1072#if INFEAS==1
1073        double sum = sumDownCost_;
1074        int number = numberTimesDown_;
1075        double originalValue = data.originalObjective_;
1076        assert (originalValue != COIN_DBL_MAX);
1077        double distanceToCutoff =  data.cutoff_  - originalValue;
1078        if (distanceToCutoff > 1.0e20)
1079            distanceToCutoff = 10.0 + fabs(originalValue);
1080        sum += numberTimesDownInfeasible_ * CoinMax(distanceToCutoff, 1.0e-12 * (1.0 + fabs(originalValue)));
1081        setDownDynamicPseudoCost(sum / static_cast<double> (number));
1082#endif
1083    } else {
1084        // up
1085        movement = ceil(value) - value;
1086        if (feasible) {
1087#ifdef COIN_DEVELOP
1088            hist.status_ = 'U';
1089#endif
1090            movement = CoinMax(movement, MINIMUM_MOVEMENT);
1091            //printf("(up change %g value down %g ",change,movement);
1092            incrementNumberTimesUp();
1093            addToSumUpChange(1.0e-30 + movement);
1094            addToSumUpDecrease(data.intDecrease_);
1095#if TYPE2==0
1096            addToSumUpCost(change / (1.0e-30 + movement));
1097            setUpDynamicPseudoCost(sumUpCost() / static_cast<double> (numberTimesUp()));
1098#elif TYPE2==1
1099            addToSumUpCost(change);
1100            setUpDynamicPseudoCost(sumUpCost() / sumUpChange());
1101#elif TYPE2==2
1102            addToSumUpCost(change*TYPERATIO + (1.0 - TYPERATIO)*change / (1.0e-30 + movement));
1103            setUpDynamicPseudoCost(sumUpCost()*(TYPERATIO / sumUpChange() + (1.0 - TYPERATIO) / (double) numberTimesUp()));
1104#endif
1105        } else {
1106#ifdef COIN_DEVELOP
1107            hist.status_ = 'u';
1108#endif
1109            //printf("(up infeasible value down %g ",change,movement);
1110            incrementNumberTimesUp();
1111            incrementNumberTimesUpInfeasible();
1112#if INFEAS==2
1113            double distanceToCutoff = 0.0;
1114            double objectiveValue = model->getCurrentMinimizationObjValue();
1115            distanceToCutoff =  model->getCutoff()  - originalValue;
1116            if (distanceToCutoff < 1.0e20)
1117                change = distanceToCutoff * 2.0;
1118            else
1119                change = upDynamicPseudoCost() * movement * 10.0;
1120            change = CoinMax(1.0e-12 * (1.0 + fabs(originalValue)), change);
1121            addToSumUpChange(1.0e-30 + movement);
1122            addToSumUpDecrease(data.intDecrease_);
1123#if TYPE2==0
1124            addToSumUpCost(change / (1.0e-30 + movement));
1125            setUpDynamicPseudoCost(sumUpCost() / (double) numberTimesUp());
1126#elif TYPE2==1
1127            addToSumUpCost(change);
1128            setUpDynamicPseudoCost(sumUpCost() / sumUpChange());
1129#elif TYPE2==2
1130            addToSumUpCost(change*TYPERATIO + (1.0 - TYPERATIO)*change / (1.0e-30 + movement));
1131            setUpDynamicPseudoCost(sumUpCost()*(TYPERATIO / sumUpChange() + (1.0 - TYPERATIO) / (double) numberTimesUp()));
1132#endif
1133#endif
1134        }
1135#if INFEAS==1
1136        double sum = sumUpCost_;
1137        int number = numberTimesUp_;
1138        double originalValue = data.originalObjective_;
1139        assert (originalValue != COIN_DBL_MAX);
1140        double distanceToCutoff =  data.cutoff_  - originalValue;
1141        if (distanceToCutoff > 1.0e20)
1142            distanceToCutoff = 10.0 + fabs(originalValue);
1143        sum += numberTimesUpInfeasible_ * CoinMax(distanceToCutoff, 1.0e-12 * (1.0 + fabs(originalValue)));
1144        setUpDynamicPseudoCost(sum / static_cast<double> (number));
1145#endif
1146    }
1147    if (data.way_ < 0)
1148        assert (numberTimesDown_ > 0);
1149    else
1150        assert (numberTimesUp_ > 0);
1151    assert (downDynamicPseudoCost_ >= 0.0 && downDynamicPseudoCost_ < 1.0e100);
1152    downDynamicPseudoCost_ = CoinMax(1.0e-10, downDynamicPseudoCost_);
1153    assert (upDynamicPseudoCost_ >= 0.0 && upDynamicPseudoCost_ < 1.0e100);
1154    upDynamicPseudoCost_ = CoinMax(1.0e-10, upDynamicPseudoCost_);
1155#ifdef COIN_DEVELOP
1156    hist.sequence_ = columnNumber_;
1157    hist.numberUp_ = numberTimesUp_;
1158    hist.numberUpInf_ = numberTimesUpInfeasible_;
1159    hist.sumUp_ = sumUpCost_;
1160    hist.upEst_ = change;
1161    hist.numberDown_ = numberTimesDown_;
1162    hist.numberDownInf_ = numberTimesDownInfeasible_;
1163    hist.sumDown_ = sumDownCost_;
1164    hist.downEst_ = movement;
1165    addRecord(hist);
1166#endif
1167    //print(1,0.5);
1168    assert (downDynamicPseudoCost_ > 1.0e-40 && upDynamicPseudoCost_ > 1.0e-40);
1169#if MOD_SHADOW>1
1170    if (upShadowPrice_ > 0.0 && numberTimesDown_ >= numberBeforeTrust_
1171            && numberTimesUp_ >= numberBeforeTrust_) {
1172        // Set negative
1173        upShadowPrice_ = -upShadowPrice_;
1174        assert (downShadowPrice_ > 0.0);
1175        downShadowPrice_ = - downShadowPrice_;
1176    }
1177#endif
1178}
1179// Updates stuff like pseudocosts after mini branch and bound
1180void
1181CbcSimpleIntegerDynamicPseudoCost::updateAfterMini(int numberDown, int numberDownInfeasible,
1182        double sumDown, int numberUp,
1183        int numberUpInfeasible, double sumUp)
1184{
1185    numberTimesDown_ = numberDown;
1186    numberTimesDownInfeasible_ = numberDownInfeasible;
1187    sumDownCost_ = sumDown;
1188    numberTimesUp_ = numberUp;
1189    numberTimesUpInfeasible_ = numberUpInfeasible;
1190    sumUpCost_ = sumUp;
1191    if (numberTimesDown_ > 0) {
1192        setDownDynamicPseudoCost(sumDownCost_ / static_cast<double> (numberTimesDown_));
1193        assert (downDynamicPseudoCost_ > 0.0 && downDynamicPseudoCost_ < 1.0e50);
1194    }
1195    if (numberTimesUp_ > 0) {
1196        setUpDynamicPseudoCost(sumUpCost_ / static_cast<double> (numberTimesUp_));
1197        assert (upDynamicPseudoCost_ > 0.0 && upDynamicPseudoCost_ < 1.0e50);
1198    }
1199    assert (downDynamicPseudoCost_ > 1.0e-40 && upDynamicPseudoCost_ > 1.0e-40);
1200}
1201// Pass in probing information
1202void
1203CbcSimpleIntegerDynamicPseudoCost::setProbingInformation(int fixedDown, int fixedUp)
1204{
1205    numberTimesProbingTotal_++;
1206    numberTimesDownLocalFixed_ = fixedDown;
1207    numberTimesDownTotalFixed_ += fixedDown;
1208    numberTimesUpLocalFixed_ = fixedUp;
1209    numberTimesUpTotalFixed_ += fixedUp;
1210}
1211// Print
1212void
1213CbcSimpleIntegerDynamicPseudoCost::print(int type, double value) const
1214{
1215    if (!type) {
1216        double meanDown = 0.0;
1217        double devDown = 0.0;
1218        if (numberTimesDown_) {
1219            meanDown = sumDownCost_ / static_cast<double> (numberTimesDown_);
1220            devDown = meanDown * meanDown - 2.0 * meanDown * sumDownCost_;
1221            if (devDown >= 0.0)
1222                devDown = sqrt(devDown);
1223        }
1224        double meanUp = 0.0;
1225        double devUp = 0.0;
1226        if (numberTimesUp_) {
1227            meanUp = sumUpCost_ / static_cast<double> (numberTimesUp_);
1228            devUp = meanUp * meanUp - 2.0 * meanUp * sumUpCost_;
1229            if (devUp >= 0.0)
1230                devUp = sqrt(devUp);
1231        }
1232        printf("%d down %d times (%d inf) mean %g (dev %g) up %d times (%d inf) mean %g (dev %g)\n",
1233               columnNumber_,
1234               numberTimesDown_, numberTimesDownInfeasible_, meanDown, devDown,
1235               numberTimesUp_, numberTimesUpInfeasible_, meanUp, devUp);
1236    } else {
1237        const double * upper = model_->getCbcColUpper();
1238        double integerTolerance =
1239            model_->getDblParam(CbcModel::CbcIntegerTolerance);
1240        double below = floor(value + integerTolerance);
1241        double above = below + 1.0;
1242        if (above > upper[columnNumber_]) {
1243            above = below;
1244            below = above - 1;
1245        }
1246        double objectiveValue = model_->getCurrentMinimizationObjValue();
1247        double distanceToCutoff =  model_->getCutoff() - objectiveValue;
1248        if (distanceToCutoff < 1.0e20)
1249            distanceToCutoff *= 10.0;
1250        else
1251            distanceToCutoff = 1.0e2 + fabs(objectiveValue);
1252        distanceToCutoff = CoinMax(distanceToCutoff, 1.0e-12 * (1.0 + fabs(objectiveValue)));
1253        double sum;
1254        int number;
1255        double downCost = CoinMax(value - below, 0.0);
1256        double downCost0 = downCost * downDynamicPseudoCost_;
1257        sum = sumDownCost();
1258        number = numberTimesDown();
1259        sum += numberTimesDownInfeasible() * (distanceToCutoff / (downCost + 1.0e-12));
1260        if (number > 0)
1261            downCost *= sum / static_cast<double> (number);
1262        else
1263            downCost  *=  downDynamicPseudoCost_;
1264        double upCost = CoinMax((above - value), 0.0);
1265        double upCost0 = upCost * upDynamicPseudoCost_;
1266        sum = sumUpCost();
1267        number = numberTimesUp();
1268        sum += numberTimesUpInfeasible() * (distanceToCutoff / (upCost + 1.0e-12));
1269        if (number > 0)
1270            upCost *= sum / static_cast<double> (number);
1271        else
1272            upCost  *=  upDynamicPseudoCost_;
1273        printf("%d down %d times %g (est %g)  up %d times %g (est %g)\n",
1274               columnNumber_,
1275               numberTimesDown_, downCost, downCost0,
1276               numberTimesUp_, upCost, upCost0);
1277    }
1278}
Note: See TracBrowser for help on using the repository browser.