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

Last change on this file since 1854 was 1854, checked in by stefan, 7 years ago

fix svn keywords property

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