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

Last change on this file since 1573 was 1573, checked in by lou, 8 years ago

Change to EPL license notice.

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