source: stable/2.8/Cbc/src/CbcSimpleIntegerDynamicPseudoCost.cpp @ 2128

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

sync with trunk rev 1901

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 51.7 KB
Line 
1// $Id: CbcSimpleIntegerDynamicPseudoCost.cpp 1902 2013-04-10 16:58:16Z forrest $
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#ifdef COIN_DEVELOP
739        char where;
740#endif
741        // was <= 10
742        //if (stateOfSearch<=1||model_->currentNode()->depth()<=-10 /* was ||maxValue>0.2*distanceToCutoff*/) {
743        if (stateOfSearch <= 2) {
744            // no branching solution
745#ifdef COIN_DEVELOP
746            where = 'i';
747#endif
748            returnValue = WEIGHT_BEFORE * minValue + (1.0 - WEIGHT_BEFORE) * maxValue;
749            if (0) {
750                double sum;
751                int number;
752                double downCost2 = CoinMax(value - below, 0.0);
753                sum = sumDownCost_;
754                number = numberTimesDown_;
755                if (number > 0)
756                    downCost2 *= sum / static_cast<double> (number);
757                else
758                    downCost2  *=  downDynamicPseudoCost_;
759                double upCost2 = CoinMax((above - value), 0.0);
760                sum = sumUpCost_;
761                number = numberTimesUp_;
762                if (number > 0)
763                    upCost2 *= sum / static_cast<double> (number);
764                else
765                    upCost2  *=  upDynamicPseudoCost_;
766                double minValue2 = CoinMin(downCost2, upCost2);
767                double maxValue2 = CoinMax(downCost2, upCost2);
768                printf("%d value %g downC %g upC %g minV %g maxV %g downC2 %g upC2 %g minV2 %g maxV2 %g\n",
769                       columnNumber_, value, downCost, upCost, minValue, maxValue,
770                       downCost2, upCost2, minValue2, maxValue2);
771            }
772        } else {
773            // some solution
774#ifdef COIN_DEVELOP
775            where = 'I';
776#endif
777#ifndef WEIGHT_PRODUCT
778            returnValue = WEIGHT_AFTER * minValue + (1.0 - WEIGHT_AFTER) * maxValue;
779#else
780            double minProductWeight = model_->getDblParam(CbcModel::CbcSmallChange);
781            returnValue = CoinMax(minValue, minProductWeight) * CoinMax(maxValue, minProductWeight);
782            //returnValue += minProductWeight*minValue;
783#endif
784        }
785        if (numberTimesUp_ < numberBeforeTrust_ ||
786                numberTimesDown_ < numberBeforeTrust_) {
787            //if (returnValue<1.0e10)
788            //returnValue += 1.0e12;
789            //else
790            returnValue *= 1.0e3;
791            if (!numberTimesUp_ && !numberTimesDown_)
792                returnValue *= 1.0e10;
793        }
794        //if (fabs(value-0.5)<1.0e-5) {
795        //returnValue = 3.0*returnValue + 0.2;
796        //} else if (value>0.9) {
797        //returnValue = 2.0*returnValue + 0.1;
798        //}
799        if (method_ == 1) {
800            // probing
801            // average
802            double up = 1.0e-15;
803            double down = 1.0e-15;
804            if (numberTimesProbingTotal_) {
805                up += numberTimesUpTotalFixed_ / static_cast<double> (numberTimesProbingTotal_);
806                down += numberTimesDownTotalFixed_ / static_cast<double> (numberTimesProbingTotal_);
807            }
808            returnValue = 1 + 10.0 * CoinMin(numberTimesDownLocalFixed_, numberTimesUpLocalFixed_) +
809                          CoinMin(down, up);
810            returnValue *= 1.0e-3;
811        }
812#ifdef COIN_DEVELOP
813        History hist;
814        hist.where_ = where;
815        hist.status_ = ' ';
816        hist.sequence_ = columnNumber_;
817        hist.numberUp_ = numberTimesUp_;
818        hist.numberUpInf_ = numberTimesUpInfeasible_;
819        hist.sumUp_ = sumUpCost_;
820        hist.upEst_ = upCost;
821        hist.numberDown_ = numberTimesDown_;
822        hist.numberDownInf_ = numberTimesDownInfeasible_;
823        hist.sumDown_ = sumDownCost_;
824        hist.downEst_ = downCost;
825        if (stateOfSearch)
826            addRecord(hist);
827#endif
828        return CoinMax(returnValue, 1.0e-15);
829    }
830}
831// Creates a branching object
832CbcBranchingObject *
833CbcSimpleIntegerDynamicPseudoCost::createCbcBranch(OsiSolverInterface * /*solver*/,
834        const OsiBranchingInformation * info, int way)
835{
836    double value = info->solution_[columnNumber_];
837    value = CoinMax(value, info->lower_[columnNumber_]);
838    value = CoinMin(value, info->upper_[columnNumber_]);
839    assert (info->upper_[columnNumber_] > info->lower_[columnNumber_]);
840    if (!info->hotstartSolution_ && priority_ != -999) {
841#ifndef NDEBUG
842        double nearest = floor(value + 0.5);
843        assert (fabs(value - nearest) > info->integerTolerance_);
844#endif
845    } else if (info->hotstartSolution_) {
846        double targetValue = info->hotstartSolution_[columnNumber_];
847        if (way > 0)
848            value = targetValue - 0.1;
849        else
850            value = targetValue + 0.1;
851    } else {
852        if (value <= info->lower_[columnNumber_])
853            value += 0.1;
854        else if (value >= info->upper_[columnNumber_])
855            value -= 0.1;
856    }
857    assert (value >= info->lower_[columnNumber_] &&
858            value <= info->upper_[columnNumber_]);
859    CbcDynamicPseudoCostBranchingObject * newObject =
860        new CbcDynamicPseudoCostBranchingObject(model_, columnNumber_, way,
861                                                value, this);
862    double up =  upDynamicPseudoCost_ * (ceil(value) - value);
863    double down =  downDynamicPseudoCost_ * (value - floor(value));
864    double changeInGuessed = up - down;
865    if (way > 0)
866        changeInGuessed = - changeInGuessed;
867    changeInGuessed = CoinMax(0.0, changeInGuessed);
868    //if (way>0)
869    //changeInGuessed += 1.0e8; // bias to stay up
870    newObject->setChangeInGuessed(changeInGuessed);
871    newObject->setOriginalObject(this);
872    return newObject;
873}
874
875// Return "up" estimate
876double
877CbcSimpleIntegerDynamicPseudoCost::upEstimate() const
878{
879    const double * solution = model_->testSolution();
880    const double * lower = model_->getCbcColLower();
881    const double * upper = model_->getCbcColUpper();
882    double value = solution[columnNumber_];
883    value = CoinMax(value, lower[columnNumber_]);
884    value = CoinMin(value, upper[columnNumber_]);
885    if (upper[columnNumber_] == lower[columnNumber_]) {
886        // fixed
887        return 0.0;
888    }
889    double integerTolerance =
890        model_->getDblParam(CbcModel::CbcIntegerTolerance);
891    double below = floor(value + integerTolerance);
892    double above = below + 1.0;
893    if (above > upper[columnNumber_]) {
894        above = below;
895        below = above - 1;
896    }
897    double upCost = CoinMax((above - value) * upDynamicPseudoCost_, 0.0);
898    return upCost;
899}
900// Return "down" estimate
901double
902CbcSimpleIntegerDynamicPseudoCost::downEstimate() const
903{
904    const double * solution = model_->testSolution();
905    const double * lower = model_->getCbcColLower();
906    const double * upper = model_->getCbcColUpper();
907    double value = solution[columnNumber_];
908    value = CoinMax(value, lower[columnNumber_]);
909    value = CoinMin(value, upper[columnNumber_]);
910    if (upper[columnNumber_] == lower[columnNumber_]) {
911        // fixed
912        return 0.0;
913    }
914    double integerTolerance =
915        model_->getDblParam(CbcModel::CbcIntegerTolerance);
916    double below = floor(value + integerTolerance);
917    double above = below + 1.0;
918    if (above > upper[columnNumber_]) {
919        above = below;
920        below = above - 1;
921    }
922    double downCost = CoinMax((value - below) * downDynamicPseudoCost_, 0.0);
923    return downCost;
924}
925// Set down pseudo cost
926void
927CbcSimpleIntegerDynamicPseudoCost::setDownDynamicPseudoCost(double value)
928{
929#ifdef TRACE_ONE
930    double oldDown = sumDownCost_;
931#endif
932    downDynamicPseudoCost_ = value;
933    sumDownCost_ = CoinMax(sumDownCost_, value * numberTimesDown_);
934#ifdef TRACE_ONE
935    if (columnNumber_ == TRACE_ONE) {
936        double down = downDynamicPseudoCost_ * numberTimesDown_;
937        printf("For %d sumDown %g (%d), inf (%d) - pseudo %g - sumDown was %g -> %g\n",
938               TRACE_ONE, down, numberTimesDown_,
939               numberTimesDownInfeasible_, downDynamicPseudoCost_,
940               oldDown, sumDownCost_);
941    }
942#endif
943}
944// Modify down pseudo cost in a slightly different way
945void
946CbcSimpleIntegerDynamicPseudoCost::updateDownDynamicPseudoCost(double value)
947{
948    sumDownCost_ += value;
949    numberTimesDown_++;
950    downDynamicPseudoCost_ = sumDownCost_ / static_cast<double>(numberTimesDown_);
951}
952// Set up pseudo cost
953void
954CbcSimpleIntegerDynamicPseudoCost::setUpDynamicPseudoCost(double value)
955{
956#ifdef TRACE_ONE
957    double oldUp = sumUpCost_;
958#endif
959    upDynamicPseudoCost_ = value;
960    sumUpCost_ = CoinMax(sumUpCost_, value * numberTimesUp_);
961#ifdef TRACE_ONE
962    if (columnNumber_ == TRACE_ONE) {
963        double up = upDynamicPseudoCost_ * numberTimesUp_;
964        printf("For %d sumUp %g (%d), inf (%d) - pseudo %g - sumUp was %g -> %g\n",
965               TRACE_ONE, up, numberTimesUp_,
966               numberTimesUpInfeasible_, upDynamicPseudoCost_,
967               oldUp, sumUpCost_);
968    }
969#endif
970}
971// Modify up pseudo cost in a slightly different way
972void
973CbcSimpleIntegerDynamicPseudoCost::updateUpDynamicPseudoCost(double value)
974{
975    sumUpCost_ += value;
976    numberTimesUp_++;
977    upDynamicPseudoCost_ = sumUpCost_ / static_cast<double>(numberTimesUp_);
978}
979/* Pass in information on branch just done and create CbcObjectUpdateData instance.
980   If object does not need data then backward pointer will be NULL.
981   Assumes can get information from solver */
982CbcObjectUpdateData
983CbcSimpleIntegerDynamicPseudoCost::createUpdateInformation(const OsiSolverInterface * solver,
984        const CbcNode * node,
985        const CbcBranchingObject * branchingObject)
986{
987    double originalValue = node->objectiveValue();
988    int originalUnsatisfied = node->numberUnsatisfied();
989    double objectiveValue = solver->getObjValue() * solver->getObjSense();
990    int unsatisfied = 0;
991    int i;
992    //might be base model - doesn't matter
993    int numberIntegers = model_->numberIntegers();;
994    const double * solution = solver->getColSolution();
995    //const double * lower = solver->getColLower();
996    //const double * upper = solver->getColUpper();
997    double change = CoinMax(0.0, objectiveValue - originalValue);
998    int iStatus;
999    if (solver->isProvenOptimal())
1000        iStatus = 0; // optimal
1001    else if (solver->isIterationLimitReached()
1002             && !solver->isDualObjectiveLimitReached())
1003        iStatus = 2; // unknown
1004    else
1005        iStatus = 1; // infeasible
1006
1007    bool feasible = iStatus != 1;
1008    if (feasible) {
1009        double integerTolerance =
1010            model_->getDblParam(CbcModel::CbcIntegerTolerance);
1011        const int * integerVariable = model_->integerVariable();
1012        for (i = 0; i < numberIntegers; i++) {
1013            int j = integerVariable[i];
1014            double value = solution[j];
1015            double nearest = floor(value + 0.5);
1016            if (fabs(value - nearest) > integerTolerance)
1017                unsatisfied++;
1018        }
1019    }
1020    int way = branchingObject->way();
1021    way = - way; // because after branch so moved on
1022    double value = branchingObject->value();
1023    CbcObjectUpdateData newData (this, way,
1024                                 change, iStatus,
1025                                 originalUnsatisfied - unsatisfied, value);
1026    newData.originalObjective_ = originalValue;
1027    // Solvers know about direction
1028    double direction = solver->getObjSense();
1029    solver->getDblParam(OsiDualObjectiveLimit, newData.cutoff_);
1030    newData.cutoff_ *= direction;
1031    return newData;
1032}
1033// Just update using feasible branches and keep count of infeasible
1034#undef INFEAS
1035// Update object by CbcObjectUpdateData
1036void
1037CbcSimpleIntegerDynamicPseudoCost::updateInformation(const CbcObjectUpdateData & data)
1038{
1039    bool feasible = data.status_ != 1;
1040    int way = data.way_;
1041    double value = data.branchingValue_;
1042    double change = data.change_;
1043#ifdef COIN_DEVELOP
1044    History hist;
1045    hist.where_ = 'U'; // need to tell if hot
1046#endif
1047    double movement = 0.0;
1048    if (way < 0) {
1049        // down
1050        movement = value - floor(value);
1051        if (feasible) {
1052#ifdef COIN_DEVELOP
1053            hist.status_ = 'D';
1054#endif
1055            movement = CoinMax(movement, MINIMUM_MOVEMENT);
1056            //printf("(down change %g value down %g ",change,movement);
1057            incrementNumberTimesDown();
1058            addToSumDownChange(1.0e-30 + movement);
1059            addToSumDownDecrease(data.intDecrease_);
1060#if TYPE2==0
1061            addToSumDownCost(change / (1.0e-30 + movement));
1062            setDownDynamicPseudoCost(sumDownCost() / static_cast<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        } else {
1071#ifdef COIN_DEVELOP
1072            hist.status_ = 'd';
1073#endif
1074            //printf("(down infeasible value down %g ",change,movement);
1075            incrementNumberTimesDown();
1076            incrementNumberTimesDownInfeasible();
1077#if INFEAS==2
1078            double distanceToCutoff = 0.0;
1079            double objectiveValue = model->getCurrentMinimizationObjValue();
1080            distanceToCutoff =  model->getCutoff()  - originalValue;
1081            if (distanceToCutoff < 1.0e20)
1082                change = distanceToCutoff * 2.0;
1083            else
1084                change = downDynamicPseudoCost() * movement * 10.0;
1085            change = CoinMax(1.0e-12 * (1.0 + fabs(originalValue)), change);
1086            addToSumDownChange(1.0e-30 + movement);
1087            addToSumDownDecrease(data.intDecrease_);
1088#if TYPE2==0
1089            addToSumDownCost(change / (1.0e-30 + movement));
1090            setDownDynamicPseudoCost(sumDownCost() / (double) numberTimesDown());
1091#elif TYPE2==1
1092            addToSumDownCost(change);
1093            setDownDynamicPseudoCost(sumDownCost() / sumDownChange());
1094#elif TYPE2==2
1095            addToSumDownCost(change*TYPERATIO + (1.0 - TYPERATIO)*change / (1.0e-30 + movement));
1096            setDownDynamicPseudoCost(sumDownCost()*(TYPERATIO / sumDownChange() + (1.0 - TYPERATIO) / (double) numberTimesDown()));
1097#endif
1098#endif
1099        }
1100#if INFEAS==1
1101        double sum = sumDownCost_;
1102        int number = numberTimesDown_;
1103        double originalValue = data.originalObjective_;
1104        assert (originalValue != COIN_DBL_MAX);
1105        double distanceToCutoff =  data.cutoff_  - originalValue;
1106        if (distanceToCutoff > 1.0e20)
1107            distanceToCutoff = 10.0 + fabs(originalValue);
1108        sum += INFEAS_MULTIPLIER*numberTimesDownInfeasible_ * CoinMax(distanceToCutoff, 1.0e-12 * (1.0 + fabs(originalValue)));
1109        setDownDynamicPseudoCost(sum / static_cast<double> (number));
1110#endif
1111    } else {
1112        // up
1113        movement = ceil(value) - value;
1114        if (feasible) {
1115#ifdef COIN_DEVELOP
1116            hist.status_ = 'U';
1117#endif
1118            movement = CoinMax(movement, MINIMUM_MOVEMENT);
1119            //printf("(up change %g value down %g ",change,movement);
1120            incrementNumberTimesUp();
1121            addToSumUpChange(1.0e-30 + movement);
1122            addToSumUpDecrease(data.intDecrease_);
1123#if TYPE2==0
1124            addToSumUpCost(change / (1.0e-30 + movement));
1125            setUpDynamicPseudoCost(sumUpCost() / static_cast<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        } else {
1134#ifdef COIN_DEVELOP
1135            hist.status_ = 'u';
1136#endif
1137            //printf("(up infeasible value down %g ",change,movement);
1138            incrementNumberTimesUp();
1139            incrementNumberTimesUpInfeasible();
1140#if INFEAS==2
1141            double distanceToCutoff = 0.0;
1142            double objectiveValue = model->getCurrentMinimizationObjValue();
1143            distanceToCutoff =  model->getCutoff()  - originalValue;
1144            if (distanceToCutoff < 1.0e20)
1145                change = distanceToCutoff * 2.0;
1146            else
1147                change = upDynamicPseudoCost() * movement * 10.0;
1148            change = CoinMax(1.0e-12 * (1.0 + fabs(originalValue)), change);
1149            addToSumUpChange(1.0e-30 + movement);
1150            addToSumUpDecrease(data.intDecrease_);
1151#if TYPE2==0
1152            addToSumUpCost(change / (1.0e-30 + movement));
1153            setUpDynamicPseudoCost(sumUpCost() / (double) numberTimesUp());
1154#elif TYPE2==1
1155            addToSumUpCost(change);
1156            setUpDynamicPseudoCost(sumUpCost() / sumUpChange());
1157#elif TYPE2==2
1158            addToSumUpCost(change*TYPERATIO + (1.0 - TYPERATIO)*change / (1.0e-30 + movement));
1159            setUpDynamicPseudoCost(sumUpCost()*(TYPERATIO / sumUpChange() + (1.0 - TYPERATIO) / (double) numberTimesUp()));
1160#endif
1161#endif
1162        }
1163#if INFEAS==1
1164        double sum = sumUpCost_;
1165        int number = numberTimesUp_;
1166        double originalValue = data.originalObjective_;
1167        assert (originalValue != COIN_DBL_MAX);
1168        double distanceToCutoff =  data.cutoff_  - originalValue;
1169        if (distanceToCutoff > 1.0e20)
1170            distanceToCutoff = 10.0 + fabs(originalValue);
1171        sum += INFEAS_MULTIPLIER*numberTimesUpInfeasible_ * CoinMax(distanceToCutoff, 1.0e-12 * (1.0 + fabs(originalValue)));
1172        setUpDynamicPseudoCost(sum / static_cast<double> (number));
1173#endif
1174    }
1175    if (data.way_ < 0)
1176        assert (numberTimesDown_ > 0);
1177    else
1178        assert (numberTimesUp_ > 0);
1179    assert (downDynamicPseudoCost_ >= 0.0 && downDynamicPseudoCost_ < 1.0e100);
1180    downDynamicPseudoCost_ = CoinMax(1.0e-10, downDynamicPseudoCost_);
1181    assert (upDynamicPseudoCost_ >= 0.0 && upDynamicPseudoCost_ < 1.0e100);
1182    upDynamicPseudoCost_ = CoinMax(1.0e-10, upDynamicPseudoCost_);
1183#ifdef COIN_DEVELOP
1184    hist.sequence_ = columnNumber_;
1185    hist.numberUp_ = numberTimesUp_;
1186    hist.numberUpInf_ = numberTimesUpInfeasible_;
1187    hist.sumUp_ = sumUpCost_;
1188    hist.upEst_ = change;
1189    hist.numberDown_ = numberTimesDown_;
1190    hist.numberDownInf_ = numberTimesDownInfeasible_;
1191    hist.sumDown_ = sumDownCost_;
1192    hist.downEst_ = movement;
1193    addRecord(hist);
1194#endif
1195    //print(1,0.5);
1196    assert (downDynamicPseudoCost_ > 1.0e-40 && upDynamicPseudoCost_ > 1.0e-40);
1197#if MOD_SHADOW>1
1198    if (upShadowPrice_ > 0.0 && numberTimesDown_ >= numberBeforeTrust_
1199            && numberTimesUp_ >= numberBeforeTrust_) {
1200        // Set negative
1201        upShadowPrice_ = -upShadowPrice_;
1202        assert (downShadowPrice_ > 0.0);
1203        downShadowPrice_ = - downShadowPrice_;
1204    }
1205#endif
1206}
1207// Updates stuff like pseudocosts after mini branch and bound
1208void
1209CbcSimpleIntegerDynamicPseudoCost::updateAfterMini(int numberDown, int numberDownInfeasible,
1210        double sumDown, int numberUp,
1211        int numberUpInfeasible, double sumUp)
1212{
1213    numberTimesDown_ = numberDown;
1214    numberTimesDownInfeasible_ = numberDownInfeasible;
1215    sumDownCost_ = sumDown;
1216    numberTimesUp_ = numberUp;
1217    numberTimesUpInfeasible_ = numberUpInfeasible;
1218    sumUpCost_ = sumUp;
1219    if (numberTimesDown_ > 0) {
1220        setDownDynamicPseudoCost(sumDownCost_ / static_cast<double> (numberTimesDown_));
1221        assert (downDynamicPseudoCost_ > 0.0 && downDynamicPseudoCost_ < 1.0e50);
1222    }
1223    if (numberTimesUp_ > 0) {
1224        setUpDynamicPseudoCost(sumUpCost_ / static_cast<double> (numberTimesUp_));
1225        assert (upDynamicPseudoCost_ > 0.0 && upDynamicPseudoCost_ < 1.0e50);
1226    }
1227    assert (downDynamicPseudoCost_ > 1.0e-40 && upDynamicPseudoCost_ > 1.0e-40);
1228}
1229// Pass in probing information
1230void
1231CbcSimpleIntegerDynamicPseudoCost::setProbingInformation(int fixedDown, int fixedUp)
1232{
1233    numberTimesProbingTotal_++;
1234    numberTimesDownLocalFixed_ = fixedDown;
1235    numberTimesDownTotalFixed_ += fixedDown;
1236    numberTimesUpLocalFixed_ = fixedUp;
1237    numberTimesUpTotalFixed_ += fixedUp;
1238}
1239// Print
1240void
1241CbcSimpleIntegerDynamicPseudoCost::print(int type, double value) const
1242{
1243    if (!type) {
1244        double meanDown = 0.0;
1245        double devDown = 0.0;
1246        if (numberTimesDown_) {
1247            meanDown = sumDownCost_ / static_cast<double> (numberTimesDown_);
1248            devDown = meanDown * meanDown - 2.0 * meanDown * sumDownCost_;
1249            if (devDown >= 0.0)
1250                devDown = sqrt(devDown);
1251        }
1252        double meanUp = 0.0;
1253        double devUp = 0.0;
1254        if (numberTimesUp_) {
1255            meanUp = sumUpCost_ / static_cast<double> (numberTimesUp_);
1256            devUp = meanUp * meanUp - 2.0 * meanUp * sumUpCost_;
1257            if (devUp >= 0.0)
1258                devUp = sqrt(devUp);
1259        }
1260        printf("%d down %d times (%d inf) mean %g (dev %g) up %d times (%d inf) mean %g (dev %g)\n",
1261               columnNumber_,
1262               numberTimesDown_, numberTimesDownInfeasible_, meanDown, devDown,
1263               numberTimesUp_, numberTimesUpInfeasible_, meanUp, devUp);
1264    } else {
1265        const double * upper = model_->getCbcColUpper();
1266        double integerTolerance =
1267            model_->getDblParam(CbcModel::CbcIntegerTolerance);
1268        double below = floor(value + integerTolerance);
1269        double above = below + 1.0;
1270        if (above > upper[columnNumber_]) {
1271            above = below;
1272            below = above - 1;
1273        }
1274        double objectiveValue = model_->getCurrentMinimizationObjValue();
1275        double distanceToCutoff =  model_->getCutoff() - objectiveValue;
1276        if (distanceToCutoff < 1.0e20)
1277            distanceToCutoff *= 10.0;
1278        else
1279            distanceToCutoff = 1.0e2 + fabs(objectiveValue);
1280        distanceToCutoff = CoinMax(distanceToCutoff, 1.0e-12 * (1.0 + fabs(objectiveValue)));
1281        double sum;
1282        int number;
1283        double downCost = CoinMax(value - below, 0.0);
1284        double downCost0 = downCost * downDynamicPseudoCost_;
1285        sum = sumDownCost();
1286        number = numberTimesDown();
1287        sum += INFEAS_MULTIPLIER*numberTimesDownInfeasible() * (distanceToCutoff / (downCost + 1.0e-12));
1288        if (number > 0)
1289            downCost *= sum / static_cast<double> (number);
1290        else
1291            downCost  *=  downDynamicPseudoCost_;
1292        double upCost = CoinMax((above - value), 0.0);
1293        double upCost0 = upCost * upDynamicPseudoCost_;
1294        sum = sumUpCost();
1295        number = numberTimesUp();
1296        sum += INFEAS_MULTIPLIER*numberTimesUpInfeasible() * (distanceToCutoff / (upCost + 1.0e-12));
1297        if (number > 0)
1298            upCost *= sum / static_cast<double> (number);
1299        else
1300            upCost  *=  upDynamicPseudoCost_;
1301        printf("%d down %d times %g (est %g)  up %d times %g (est %g)\n",
1302               columnNumber_,
1303               numberTimesDown_, downCost, downCost0,
1304               numberTimesUp_, upCost, upCost0);
1305    }
1306}
1307
1308//##############################################################################
1309
1310// Default Constructor
1311CbcIntegerPseudoCostBranchingObject::CbcIntegerPseudoCostBranchingObject()
1312        : CbcIntegerBranchingObject()
1313{
1314    changeInGuessed_ = 1.0e-5;
1315}
1316
1317// Useful constructor
1318CbcIntegerPseudoCostBranchingObject::CbcIntegerPseudoCostBranchingObject (CbcModel * model,
1319        int variable, int way , double value)
1320        : CbcIntegerBranchingObject(model, variable, way, value)
1321{
1322}
1323// Useful constructor for fixing
1324CbcIntegerPseudoCostBranchingObject::CbcIntegerPseudoCostBranchingObject (CbcModel * model,
1325        int variable, int way,
1326        double lowerValue,
1327        double /*upperValue*/)
1328        : CbcIntegerBranchingObject(model, variable, way, lowerValue)
1329{
1330    changeInGuessed_ = 1.0e100;
1331}
1332
1333
1334// Copy constructor
1335CbcIntegerPseudoCostBranchingObject::CbcIntegerPseudoCostBranchingObject (
1336    const CbcIntegerPseudoCostBranchingObject & rhs)
1337        : CbcIntegerBranchingObject(rhs)
1338{
1339    changeInGuessed_ = rhs.changeInGuessed_;
1340}
1341
1342// Assignment operator
1343CbcIntegerPseudoCostBranchingObject &
1344CbcIntegerPseudoCostBranchingObject::operator=( const CbcIntegerPseudoCostBranchingObject & rhs)
1345{
1346    if (this != &rhs) {
1347        CbcIntegerBranchingObject::operator=(rhs);
1348        changeInGuessed_ = rhs.changeInGuessed_;
1349    }
1350    return *this;
1351}
1352CbcBranchingObject *
1353CbcIntegerPseudoCostBranchingObject::clone() const
1354{
1355    return (new CbcIntegerPseudoCostBranchingObject(*this));
1356}
1357
1358
1359// Destructor
1360CbcIntegerPseudoCostBranchingObject::~CbcIntegerPseudoCostBranchingObject ()
1361{
1362}
1363
1364/*
1365  Perform a branch by adjusting the bounds of the specified variable. Note
1366  that each arm of the branch advances the object to the next arm by
1367  advancing the value of way_.
1368
1369  Providing new values for the variable's lower and upper bounds for each
1370  branching direction gives a little bit of additional flexibility and will
1371  be easily extensible to multi-way branching.
1372  Returns change in guessed objective on next branch
1373*/
1374double
1375CbcIntegerPseudoCostBranchingObject::branch()
1376{
1377    CbcIntegerBranchingObject::branch();
1378    return changeInGuessed_;
1379}
1380
1381/** Compare the \c this with \c brObj. \c this and \c brObj must be os the
1382    same type and must have the same original object, but they may have
1383    different feasible regions.
1384    Return the appropriate CbcRangeCompare value (first argument being the
1385    sub/superset if that's the case). In case of overlap (and if \c
1386    replaceIfOverlap is true) replace the current branching object with one
1387    whose feasible region is the overlap.
1388*/
1389CbcRangeCompare
1390CbcIntegerPseudoCostBranchingObject::compareBranchingObject
1391(const CbcBranchingObject* brObj, const bool replaceIfOverlap)
1392{
1393    const CbcIntegerPseudoCostBranchingObject* br =
1394        dynamic_cast<const CbcIntegerPseudoCostBranchingObject*>(brObj);
1395    assert(br);
1396    double* thisBd = way_ < 0 ? down_ : up_;
1397    const double* otherBd = br->way_ < 0 ? br->down_ : br->up_;
1398    return CbcCompareRanges(thisBd, otherBd, replaceIfOverlap);
1399}
1400
Note: See TracBrowser for help on using the repository browser.