source: stable/2.9/Cbc/src/CbcSimpleIntegerDynamicPseudoCost.cpp @ 2180

Last change on this file since 2180 was 1943, checked in by forrest, 6 years ago

more options, copy statistics structure analysis
start coding of "switch" variables i.e. badly scaled ints or hi/lo
changes to allow more influence on small branch and bound
changes to get correct printout with threads
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 65.5 KB
Line 
1// $Id: CbcSimpleIntegerDynamicPseudoCost.cpp 1943 2013-07-21 09:05:45Z 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_HAS_CLP
27#include "OsiClpSolverInterface.hpp"
28#endif
29#ifdef COIN_DEVELOP
30typedef struct {
31    double sumUp_;
32    double upEst_; // or change in obj in update
33    double sumDown_;
34    double downEst_; // or movement in value in update
35    int sequence_;
36    int numberUp_;
37    int numberUpInf_;
38    int numberDown_;
39    int numberDownInf_;
40    char where_;
41    char status_;
42} History;
43History * history = NULL;
44int numberHistory = 0;
45int maxHistory = 0;
46bool getHistoryStatistics_ = true;
47static void increaseHistory()
48{
49    if (numberHistory == maxHistory) {
50        maxHistory = 100 + (3 * maxHistory) / 2;
51        History * temp = new History [maxHistory];
52        memcpy(temp, history, numberHistory*sizeof(History));
53        delete [] history;
54        history = temp;
55    }
56}
57static bool addRecord(History newOne)
58{
59    //if (!getHistoryStatistics_)
60    return false;
61    bool fromCompare = false;
62    int i;
63    for ( i = numberHistory - 1; i >= 0; i--) {
64        if (newOne.sequence_ != history[i].sequence_)
65            continue;
66        if (newOne.where_ != history[i].where_)
67            continue;
68        if (newOne.numberUp_ != history[i].numberUp_)
69            continue;
70        if (newOne.sumUp_ != history[i].sumUp_)
71            continue;
72        if (newOne.numberUpInf_ != history[i].numberUpInf_)
73            continue;
74        if (newOne.upEst_ != history[i].upEst_)
75            continue;
76        if (newOne.numberDown_ != history[i].numberDown_)
77            continue;
78        if (newOne.sumDown_ != history[i].sumDown_)
79            continue;
80        if (newOne.numberDownInf_ != history[i].numberDownInf_)
81            continue;
82        if (newOne.downEst_ != history[i].downEst_)
83            continue;
84        // If B knock out previous B
85        if (newOne.where_ == 'C') {
86            fromCompare = true;
87            if (newOne.status_ == 'B') {
88                int j;
89                for (j = i - 1; j >= 0; j--) {
90                    if (history[j].where_ == 'C') {
91                        if (history[j].status_ == 'I') {
92                            break;
93                        } else if (history[j].status_ == 'B') {
94                            history[j].status_ = ' ';
95                            break;
96                        }
97                    }
98                }
99            }
100            break;
101        }
102    }
103    if (i == -1 || fromCompare) {
104        //add
105        increaseHistory();
106        history[numberHistory++] = newOne;
107        return true;
108    } else {
109        return false;
110    }
111}
112#endif
113
114/** Default Constructor
115
116  Equivalent to an unspecified binary variable.
117*/
118CbcSimpleIntegerDynamicPseudoCost::CbcSimpleIntegerDynamicPseudoCost ()
119        : CbcSimpleInteger(),
120        downDynamicPseudoCost_(1.0e-5),
121        upDynamicPseudoCost_(1.0e-5),
122        upDownSeparator_(-1.0),
123        sumDownCost_(0.0),
124        sumUpCost_(0.0),
125        sumDownChange_(0.0),
126        sumUpChange_(0.0),
127        downShadowPrice_(0.0),
128        upShadowPrice_(0.0),
129        sumDownDecrease_(0.0),
130        sumUpDecrease_(0.0),
131        lastDownCost_(0.0),
132        lastUpCost_(0.0),
133        lastDownDecrease_(0),
134        lastUpDecrease_(0),
135        numberTimesDown_(0),
136        numberTimesUp_(0),
137        numberTimesDownInfeasible_(0),
138        numberTimesUpInfeasible_(0),
139        numberBeforeTrust_(0),
140        numberTimesDownLocalFixed_(0),
141        numberTimesUpLocalFixed_(0),
142        numberTimesDownTotalFixed_(0.0),
143        numberTimesUpTotalFixed_(0.0),
144        numberTimesProbingTotal_(0),
145        method_(0)
146{
147}
148
149/** Useful constructor
150
151  Loads dynamic upper & lower bounds for the specified variable.
152*/
153CbcSimpleIntegerDynamicPseudoCost::CbcSimpleIntegerDynamicPseudoCost (CbcModel * model,
154        int iColumn, double breakEven)
155        : CbcSimpleInteger(model, iColumn, breakEven),
156        upDownSeparator_(-1.0),
157        sumDownCost_(0.0),
158        sumUpCost_(0.0),
159        sumDownChange_(0.0),
160        sumUpChange_(0.0),
161        downShadowPrice_(0.0),
162        upShadowPrice_(0.0),
163        sumDownDecrease_(0.0),
164        sumUpDecrease_(0.0),
165        lastDownCost_(0.0),
166        lastUpCost_(0.0),
167        lastDownDecrease_(0),
168        lastUpDecrease_(0),
169        numberTimesDown_(0),
170        numberTimesUp_(0),
171        numberTimesDownInfeasible_(0),
172        numberTimesUpInfeasible_(0),
173        numberBeforeTrust_(0),
174        numberTimesDownLocalFixed_(0),
175        numberTimesUpLocalFixed_(0),
176        numberTimesDownTotalFixed_(0.0),
177        numberTimesUpTotalFixed_(0.0),
178        numberTimesProbingTotal_(0),
179        method_(0)
180{
181    const double * cost = model->getObjCoefficients();
182    double costValue = CoinMax(1.0e-5, fabs(cost[iColumn]));
183    // treat as if will cost what it says up
184    upDynamicPseudoCost_ = costValue;
185    // and balance at breakeven
186    downDynamicPseudoCost_ = ((1.0 - breakEven_) * upDynamicPseudoCost_) / breakEven_;
187    // so initial will have some effect
188    sumUpCost_ = 2.0 * upDynamicPseudoCost_;
189    sumUpChange_ = 2.0;
190    numberTimesUp_ = 2;
191    sumDownCost_ = 2.0 * downDynamicPseudoCost_;
192    sumDownChange_ = 2.0;
193    numberTimesDown_ = 2;
194#if TYPE2==0
195    // No
196    sumUpCost_ = 0.0;
197    sumUpChange_ = 0.0;
198    numberTimesUp_ = 0;
199    sumDownCost_ = 0.0;
200    sumDownChange_ = 0.0;
201    numberTimesDown_ = 0;
202#else
203    sumUpCost_ = 1.0 * upDynamicPseudoCost_;
204    sumUpChange_ = 1.0;
205    numberTimesUp_ = 1;
206    sumDownCost_ = 1.0 * downDynamicPseudoCost_;
207    sumDownChange_ = 1.0;
208    numberTimesDown_ = 1;
209#endif
210}
211
212/** Useful constructor
213
214  Loads dynamic upper & lower bounds for the specified variable.
215*/
216CbcSimpleIntegerDynamicPseudoCost::CbcSimpleIntegerDynamicPseudoCost (CbcModel * model,
217        int iColumn, double downDynamicPseudoCost,
218        double upDynamicPseudoCost)
219        : CbcSimpleInteger(model, iColumn),
220        upDownSeparator_(-1.0),
221        sumDownCost_(0.0),
222        sumUpCost_(0.0),
223        sumDownChange_(0.0),
224        sumUpChange_(0.0),
225        downShadowPrice_(0.0),
226        upShadowPrice_(0.0),
227        sumDownDecrease_(0.0),
228        sumUpDecrease_(0.0),
229        lastDownCost_(0.0),
230        lastUpCost_(0.0),
231        lastDownDecrease_(0),
232        lastUpDecrease_(0),
233        numberTimesDown_(0),
234        numberTimesUp_(0),
235        numberTimesDownInfeasible_(0),
236        numberTimesUpInfeasible_(0),
237        numberBeforeTrust_(0),
238        numberTimesDownLocalFixed_(0),
239        numberTimesUpLocalFixed_(0),
240        numberTimesDownTotalFixed_(0.0),
241        numberTimesUpTotalFixed_(0.0),
242        numberTimesProbingTotal_(0),
243        method_(0)
244{
245    downDynamicPseudoCost_ = downDynamicPseudoCost;
246    upDynamicPseudoCost_ = upDynamicPseudoCost;
247    breakEven_ = upDynamicPseudoCost_ / (upDynamicPseudoCost_ + downDynamicPseudoCost_);
248    // so initial will have some effect
249    sumUpCost_ = 2.0 * upDynamicPseudoCost_;
250    sumUpChange_ = 2.0;
251    numberTimesUp_ = 2;
252    sumDownCost_ = 2.0 * downDynamicPseudoCost_;
253    sumDownChange_ = 2.0;
254    numberTimesDown_ = 2;
255#if TYPE2==0
256    // No
257    sumUpCost_ = 0.0;
258    sumUpChange_ = 0.0;
259    numberTimesUp_ = 0;
260    sumDownCost_ = 0.0;
261    sumDownChange_ = 0.0;
262    numberTimesDown_ = 0;
263    sumUpCost_ = 1.0e-4 * upDynamicPseudoCost_;
264    sumDownCost_ = 1.0e-4 * downDynamicPseudoCost_;
265#else
266    sumUpCost_ = 1.0 * upDynamicPseudoCost_;
267    sumUpChange_ = 1.0;
268    numberTimesUp_ = 1;
269    sumDownCost_ = 1.0 * downDynamicPseudoCost_;
270    sumDownChange_ = 1.0;
271    numberTimesDown_ = 1;
272#endif
273}
274/** Useful constructor
275
276  Loads dynamic upper & lower bounds for the specified variable.
277*/
278CbcSimpleIntegerDynamicPseudoCost::CbcSimpleIntegerDynamicPseudoCost (CbcModel * model,
279        int /*dummy*/,
280        int iColumn, double downDynamicPseudoCost,
281        double upDynamicPseudoCost)
282{
283    CbcSimpleIntegerDynamicPseudoCost(model, iColumn, downDynamicPseudoCost, upDynamicPseudoCost);
284}
285
286// Copy constructor
287CbcSimpleIntegerDynamicPseudoCost::CbcSimpleIntegerDynamicPseudoCost ( const CbcSimpleIntegerDynamicPseudoCost & rhs)
288        : CbcSimpleInteger(rhs),
289        downDynamicPseudoCost_(rhs.downDynamicPseudoCost_),
290        upDynamicPseudoCost_(rhs.upDynamicPseudoCost_),
291        upDownSeparator_(rhs.upDownSeparator_),
292        sumDownCost_(rhs.sumDownCost_),
293        sumUpCost_(rhs.sumUpCost_),
294        sumDownChange_(rhs.sumDownChange_),
295        sumUpChange_(rhs.sumUpChange_),
296        downShadowPrice_(rhs.downShadowPrice_),
297        upShadowPrice_(rhs.upShadowPrice_),
298        sumDownDecrease_(rhs.sumDownDecrease_),
299        sumUpDecrease_(rhs.sumUpDecrease_),
300        lastDownCost_(rhs.lastDownCost_),
301        lastUpCost_(rhs.lastUpCost_),
302        lastDownDecrease_(rhs.lastDownDecrease_),
303        lastUpDecrease_(rhs.lastUpDecrease_),
304        numberTimesDown_(rhs.numberTimesDown_),
305        numberTimesUp_(rhs.numberTimesUp_),
306        numberTimesDownInfeasible_(rhs.numberTimesDownInfeasible_),
307        numberTimesUpInfeasible_(rhs.numberTimesUpInfeasible_),
308        numberBeforeTrust_(rhs.numberBeforeTrust_),
309        numberTimesDownLocalFixed_(rhs.numberTimesDownLocalFixed_),
310        numberTimesUpLocalFixed_(rhs.numberTimesUpLocalFixed_),
311        numberTimesDownTotalFixed_(rhs.numberTimesDownTotalFixed_),
312        numberTimesUpTotalFixed_(rhs.numberTimesUpTotalFixed_),
313        numberTimesProbingTotal_(rhs.numberTimesProbingTotal_),
314        method_(rhs.method_)
315
316{
317}
318
319// Clone
320CbcObject *
321CbcSimpleIntegerDynamicPseudoCost::clone() const
322{
323    return new CbcSimpleIntegerDynamicPseudoCost(*this);
324}
325
326// Assignment operator
327CbcSimpleIntegerDynamicPseudoCost &
328CbcSimpleIntegerDynamicPseudoCost::operator=( const CbcSimpleIntegerDynamicPseudoCost & rhs)
329{
330    if (this != &rhs) {
331        CbcSimpleInteger::operator=(rhs);
332        downDynamicPseudoCost_ = rhs.downDynamicPseudoCost_;
333        upDynamicPseudoCost_ = rhs.upDynamicPseudoCost_;
334        upDownSeparator_ = rhs.upDownSeparator_;
335        sumDownCost_ = rhs.sumDownCost_;
336        sumUpCost_ = rhs.sumUpCost_;
337        sumDownChange_ = rhs.sumDownChange_;
338        sumUpChange_ = rhs.sumUpChange_;
339        downShadowPrice_ = rhs.downShadowPrice_;
340        upShadowPrice_ = rhs.upShadowPrice_;
341        sumDownDecrease_ = rhs.sumDownDecrease_;
342        sumUpDecrease_ = rhs.sumUpDecrease_;
343        lastDownCost_ = rhs.lastDownCost_;
344        lastUpCost_ = rhs.lastUpCost_;
345        lastDownDecrease_ = rhs.lastDownDecrease_;
346        lastUpDecrease_ = rhs.lastUpDecrease_;
347        numberTimesDown_ = rhs.numberTimesDown_;
348        numberTimesUp_ = rhs.numberTimesUp_;
349        numberTimesDownInfeasible_ = rhs.numberTimesDownInfeasible_;
350        numberTimesUpInfeasible_ = rhs.numberTimesUpInfeasible_;
351        numberBeforeTrust_ = rhs.numberBeforeTrust_;
352        numberTimesDownLocalFixed_ = rhs.numberTimesDownLocalFixed_;
353        numberTimesUpLocalFixed_ = rhs.numberTimesUpLocalFixed_;
354        numberTimesDownTotalFixed_ = rhs.numberTimesDownTotalFixed_;
355        numberTimesUpTotalFixed_ = rhs.numberTimesUpTotalFixed_;
356        numberTimesProbingTotal_ = rhs.numberTimesProbingTotal_;
357        method_ = rhs.method_;
358    }
359    return *this;
360}
361
362// Destructor
363CbcSimpleIntegerDynamicPseudoCost::~CbcSimpleIntegerDynamicPseudoCost ()
364{
365}
366// Copy some information i.e. just variable stuff
367void
368CbcSimpleIntegerDynamicPseudoCost::copySome(const CbcSimpleIntegerDynamicPseudoCost * otherObject)
369{
370    downDynamicPseudoCost_ = otherObject->downDynamicPseudoCost_;
371    upDynamicPseudoCost_ = otherObject->upDynamicPseudoCost_;
372    sumDownCost_ = otherObject->sumDownCost_;
373    sumUpCost_ = otherObject->sumUpCost_;
374    sumDownChange_ = otherObject->sumDownChange_;
375    sumUpChange_ = otherObject->sumUpChange_;
376    downShadowPrice_ = otherObject->downShadowPrice_;
377    upShadowPrice_ = otherObject->upShadowPrice_;
378    sumDownDecrease_ = otherObject->sumDownDecrease_;
379    sumUpDecrease_ = otherObject->sumUpDecrease_;
380    lastDownCost_ = otherObject->lastDownCost_;
381    lastUpCost_ = otherObject->lastUpCost_;
382    lastDownDecrease_ = otherObject->lastDownDecrease_;
383    lastUpDecrease_ = otherObject->lastUpDecrease_;
384    numberTimesDown_ = otherObject->numberTimesDown_;
385    numberTimesUp_ = otherObject->numberTimesUp_;
386    numberTimesDownInfeasible_ = otherObject->numberTimesDownInfeasible_;
387    numberTimesUpInfeasible_ = otherObject->numberTimesUpInfeasible_;
388    numberTimesDownLocalFixed_ = otherObject->numberTimesDownLocalFixed_;
389    numberTimesUpLocalFixed_ = otherObject->numberTimesUpLocalFixed_;
390    numberTimesDownTotalFixed_ = otherObject->numberTimesDownTotalFixed_;
391    numberTimesUpTotalFixed_ = otherObject->numberTimesUpTotalFixed_;
392    numberTimesProbingTotal_ = otherObject->numberTimesProbingTotal_;
393}
394// Updates stuff like pseudocosts before threads
395void
396CbcSimpleIntegerDynamicPseudoCost::updateBefore(const OsiObject * rhs)
397{
398#ifndef NDEBUG
399    const CbcSimpleIntegerDynamicPseudoCost * rhsObject =
400        dynamic_cast <const CbcSimpleIntegerDynamicPseudoCost *>(rhs) ;
401    assert (rhsObject);
402#else
403    const CbcSimpleIntegerDynamicPseudoCost * rhsObject =
404        static_cast <const CbcSimpleIntegerDynamicPseudoCost *>(rhs) ;
405#endif
406    copySome(rhsObject);
407}
408// Updates stuff like pseudocosts after threads finished
409void
410CbcSimpleIntegerDynamicPseudoCost::updateAfter(const OsiObject * rhs, const OsiObject * baseObjectX)
411{
412#ifndef NDEBUG
413    const CbcSimpleIntegerDynamicPseudoCost * rhsObject =
414        dynamic_cast <const CbcSimpleIntegerDynamicPseudoCost *>(rhs) ;
415    assert (rhsObject);
416    const CbcSimpleIntegerDynamicPseudoCost * baseObject =
417        dynamic_cast <const CbcSimpleIntegerDynamicPseudoCost *>(baseObjectX) ;
418    assert (baseObject);
419#else
420    const CbcSimpleIntegerDynamicPseudoCost * rhsObject =
421        static_cast <const CbcSimpleIntegerDynamicPseudoCost *>(rhs) ;
422    const CbcSimpleIntegerDynamicPseudoCost * baseObject =
423        static_cast <const CbcSimpleIntegerDynamicPseudoCost *>(baseObjectX) ;
424#endif
425    // compute current
426    double sumDown = downDynamicPseudoCost_ * numberTimesDown_;
427    sumDown -= baseObject->downDynamicPseudoCost_ * baseObject->numberTimesDown_;
428    sumDown = CoinMax(sumDown, 0.0);
429    sumDown += rhsObject->downDynamicPseudoCost_ * rhsObject->numberTimesDown_;
430    assert (rhsObject->numberTimesDown_ >= baseObject->numberTimesDown_);
431    assert (rhsObject->numberTimesDownInfeasible_ >= baseObject->numberTimesDownInfeasible_);
432    assert( rhsObject->sumDownCost_ >= baseObject->sumDownCost_);
433    double sumUp = upDynamicPseudoCost_ * numberTimesUp_;
434    sumUp -= baseObject->upDynamicPseudoCost_ * baseObject->numberTimesUp_;
435    sumUp = CoinMax(sumUp, 0.0);
436    sumUp += rhsObject->upDynamicPseudoCost_ * rhsObject->numberTimesUp_;
437    assert (rhsObject->numberTimesUp_ >= baseObject->numberTimesUp_);
438    assert (rhsObject->numberTimesUpInfeasible_ >= baseObject->numberTimesUpInfeasible_);
439    assert( rhsObject->sumUpCost_ >= baseObject->sumUpCost_);
440    sumDownCost_ += rhsObject->sumDownCost_ - baseObject->sumDownCost_;
441    sumUpCost_ += rhsObject->sumUpCost_ - baseObject->sumUpCost_;
442    sumDownChange_ += rhsObject->sumDownChange_ - baseObject->sumDownChange_;
443    sumUpChange_ += rhsObject->sumUpChange_ - baseObject->sumUpChange_;
444    downShadowPrice_ = 0.0;
445    upShadowPrice_ = 0.0;
446    sumDownDecrease_ += rhsObject->sumDownDecrease_ - baseObject->sumDownDecrease_;
447    sumUpDecrease_ += rhsObject->sumUpDecrease_ - baseObject->sumUpDecrease_;
448    lastDownCost_ += rhsObject->lastDownCost_ - baseObject->lastDownCost_;
449    lastUpCost_ += rhsObject->lastUpCost_ - baseObject->lastUpCost_;
450    lastDownDecrease_ += rhsObject->lastDownDecrease_ - baseObject->lastDownDecrease_;
451    lastUpDecrease_ += rhsObject->lastUpDecrease_ - baseObject->lastUpDecrease_;
452    numberTimesDown_ += rhsObject->numberTimesDown_ - baseObject->numberTimesDown_;
453    numberTimesUp_ += rhsObject->numberTimesUp_ - baseObject->numberTimesUp_;
454    numberTimesDownInfeasible_ += rhsObject->numberTimesDownInfeasible_ - baseObject->numberTimesDownInfeasible_;
455    numberTimesUpInfeasible_ += rhsObject->numberTimesUpInfeasible_ - baseObject->numberTimesUpInfeasible_;
456    numberTimesDownLocalFixed_ += rhsObject->numberTimesDownLocalFixed_ - baseObject->numberTimesDownLocalFixed_;
457    numberTimesUpLocalFixed_ += rhsObject->numberTimesUpLocalFixed_ - baseObject->numberTimesUpLocalFixed_;
458    numberTimesDownTotalFixed_ += rhsObject->numberTimesDownTotalFixed_ - baseObject->numberTimesDownTotalFixed_;
459    numberTimesUpTotalFixed_ += rhsObject->numberTimesUpTotalFixed_ - baseObject->numberTimesUpTotalFixed_;
460    numberTimesProbingTotal_ += rhsObject->numberTimesProbingTotal_ - baseObject->numberTimesProbingTotal_;
461    if (numberTimesDown_ > 0) {
462        setDownDynamicPseudoCost(sumDown / static_cast<double> (numberTimesDown_));
463    }
464    if (numberTimesUp_ > 0) {
465        setUpDynamicPseudoCost(sumUp / static_cast<double> (numberTimesUp_));
466    }
467    //printf("XX %d down %d %d %g up %d %d %g\n",columnNumber_,numberTimesDown_,numberTimesDownInfeasible_,downDynamicPseudoCost_,
468    // numberTimesUp_,numberTimesUpInfeasible_,upDynamicPseudoCost_);
469    assert (downDynamicPseudoCost_ > 1.0e-40 && upDynamicPseudoCost_ > 1.0e-40);
470}
471// Same - returns true if contents match(ish)
472bool
473CbcSimpleIntegerDynamicPseudoCost::same(const CbcSimpleIntegerDynamicPseudoCost * otherObject) const
474{
475    bool okay = true;
476    if (downDynamicPseudoCost_ != otherObject->downDynamicPseudoCost_)
477        okay = false;
478    if (upDynamicPseudoCost_ != otherObject->upDynamicPseudoCost_)
479        okay = false;
480    if (sumDownCost_ != otherObject->sumDownCost_)
481        okay = false;
482    if (sumUpCost_ != otherObject->sumUpCost_)
483        okay = false;
484    if (sumDownChange_ != otherObject->sumDownChange_)
485        okay = false;
486    if (sumUpChange_ != otherObject->sumUpChange_)
487        okay = false;
488    if (downShadowPrice_ != otherObject->downShadowPrice_)
489        okay = false;
490    if (upShadowPrice_ != otherObject->upShadowPrice_)
491        okay = false;
492    if (sumDownDecrease_ != otherObject->sumDownDecrease_)
493        okay = false;
494    if (sumUpDecrease_ != otherObject->sumUpDecrease_)
495        okay = false;
496    if (lastDownCost_ != otherObject->lastDownCost_)
497        okay = false;
498    if (lastUpCost_ != otherObject->lastUpCost_)
499        okay = false;
500    if (lastDownDecrease_ != otherObject->lastDownDecrease_)
501        okay = false;
502    if (lastUpDecrease_ != otherObject->lastUpDecrease_)
503        okay = false;
504    if (numberTimesDown_ != otherObject->numberTimesDown_)
505        okay = false;
506    if (numberTimesUp_ != otherObject->numberTimesUp_)
507        okay = false;
508    if (numberTimesDownInfeasible_ != otherObject->numberTimesDownInfeasible_)
509        okay = false;
510    if (numberTimesUpInfeasible_ != otherObject->numberTimesUpInfeasible_)
511        okay = false;
512    if (numberTimesDownLocalFixed_ != otherObject->numberTimesDownLocalFixed_)
513        okay = false;
514    if (numberTimesUpLocalFixed_ != otherObject->numberTimesUpLocalFixed_)
515        okay = false;
516    if (numberTimesDownTotalFixed_ != otherObject->numberTimesDownTotalFixed_)
517        okay = false;
518    if (numberTimesUpTotalFixed_ != otherObject->numberTimesUpTotalFixed_)
519        okay = false;
520    if (numberTimesProbingTotal_ != otherObject->numberTimesProbingTotal_)
521        okay = false;
522    return okay;
523}
524/* Create an OsiSolverBranch object
525
526This returns NULL if branch not represented by bound changes
527*/
528OsiSolverBranch *
529CbcSimpleIntegerDynamicPseudoCost::solverBranch() const
530{
531    OsiSolverInterface * solver = model_->solver();
532    const double * solution = model_->testSolution();
533    const double * lower = solver->getColLower();
534    const double * upper = solver->getColUpper();
535    double value = solution[columnNumber_];
536    value = CoinMax(value, lower[columnNumber_]);
537    value = CoinMin(value, upper[columnNumber_]);
538    assert (upper[columnNumber_] > lower[columnNumber_]);
539#ifndef NDEBUG
540    double nearest = floor(value + 0.5);
541    double integerTolerance =
542        model_->getDblParam(CbcModel::CbcIntegerTolerance);
543    assert (fabs(value - nearest) > integerTolerance);
544#endif
545    OsiSolverBranch * branch = new OsiSolverBranch();
546    branch->addBranch(columnNumber_, value);
547    return branch;
548}
549//#define FUNNY_BRANCHING
550double
551CbcSimpleIntegerDynamicPseudoCost::infeasibility(const OsiBranchingInformation * info,
552        int &preferredWay) const
553{
554    assert (downDynamicPseudoCost_ > 1.0e-40 && upDynamicPseudoCost_ > 1.0e-40);
555    const double * solution = model_->testSolution();
556    const double * lower = model_->getCbcColLower();
557    const double * upper = model_->getCbcColUpper();
558#ifdef FUNNY_BRANCHING2
559    const double * dj = model_->getCbcReducedCost();
560    double djValue = dj[columnNumber_];
561    lastDownDecrease_++;
562    if (djValue > 1.0e-6) {
563        // wants to go down
564        if (true || lower[columnNumber_] > originalLower_) {
565            // Lower bound active
566            lastUpDecrease_++;
567        }
568    } else if (djValue < -1.0e-6) {
569        // wants to go up
570        if (true || upper[columnNumber_] < originalUpper_) {
571            // Upper bound active
572            lastUpDecrease_++;
573        }
574    }
575#endif
576    if (upper[columnNumber_] == lower[columnNumber_]) {
577        // fixed
578        preferredWay = 1;
579        return 0.0;
580    }
581    assert (breakEven_ > 0.0 && breakEven_ < 1.0);
582        /*
583          Find nearest integer, and integers above and below current value.
584
585          Given that we've already forced value within bounds, if
586          (current value)+(integer tolerance) > (upper bound)
587          shouldn't we declare this variable integer?
588        */
589
590    double value = solution[columnNumber_];
591    value = CoinMax(value, lower[columnNumber_]);
592    value = CoinMin(value, upper[columnNumber_]);
593    /*printf("%d %g %g %g %g\n",columnNumber_,value,lower[columnNumber_],
594      solution[columnNumber_],upper[columnNumber_]);*/
595    double nearest = floor(value + 0.5);
596    double integerTolerance =
597        model_->getDblParam(CbcModel::CbcIntegerTolerance);
598    double below = floor(value + integerTolerance);
599    double above = below + 1.0;
600    if (above > upper[columnNumber_]) {
601        above = below;
602        below = above - 1;
603    }
604#if INFEAS==1
605/*
606  Why do we inflate the distance to the cutoff by a factor of 10 for
607  values that could be considered reachable? Why do we add 100 for values
608  larger than 1e20?
609*/
610    double distanceToCutoff = 0.0;
611    double objectiveValue = model_->getCurrentMinimizationObjValue();
612    distanceToCutoff =  model_->getCutoff()  - objectiveValue;
613    if (distanceToCutoff < 1.0e20)
614        distanceToCutoff *= 10.0;
615    else
616        distanceToCutoff = 1.0e2 + fabs(objectiveValue);
617    distanceToCutoff = CoinMax(distanceToCutoff, 1.0e-12 * (1.0 + fabs(objectiveValue)));
618#endif
619    double sum;
620#ifndef INFEAS_MULTIPLIER
621#define INFEAS_MULTIPLIER 1.0
622#endif
623    double number;
624    double downCost = CoinMax(value - below, 0.0);
625#if TYPE2==0
626    sum = sumDownCost_;
627    number = numberTimesDown_;
628#if INFEAS==1
629    sum += INFEAS_MULTIPLIER*numberTimesDownInfeasible_ * CoinMax(distanceToCutoff / (downCost + 1.0e-12), sumDownCost_);
630#endif
631#elif TYPE2==1
632    sum = sumDownCost_;
633    number = sumDownChange_;
634#if INFEAS==1
635    sum += INFEAS_MULTIPLIER*numberTimesDownInfeasible_ * CoinMax(distanceToCutoff / (downCost + 1.0e-12), sumDownCost_);
636#endif
637#elif TYPE2==2
638    abort();
639#if INFEAS==1
640    sum += INFEAS_MULTIPLIER*numberTimesDownInfeasible_ * (distanceToCutoff / (downCost + 1.0e-12));
641#endif
642#endif
643#if MOD_SHADOW>0
644    if (!downShadowPrice_) {
645        if (number > 0.0)
646            downCost *= sum / number;
647        else
648            downCost  *=  downDynamicPseudoCost_;
649    } else if (downShadowPrice_ > 0.0) {
650        downCost *= downShadowPrice_;
651    } else {
652        downCost *= (downDynamicPseudoCost_ - downShadowPrice_);
653    }
654#else
655    if (downShadowPrice_ <= 0.0) {
656        if (number > 0.0)
657            downCost *= sum / number;
658        else
659            downCost  *=  downDynamicPseudoCost_;
660    } else {
661        downCost *= downShadowPrice_;
662    }
663#endif
664    double upCost = CoinMax((above - value), 0.0);
665#if TYPE2==0
666    sum = sumUpCost_;
667    number = numberTimesUp_;
668#if INFEAS==1
669    sum += INFEAS_MULTIPLIER*numberTimesUpInfeasible_ * CoinMax(distanceToCutoff / (upCost + 1.0e-12), sumUpCost_);
670#endif
671#elif TYPE2==1
672    sum = sumUpCost_;
673    number = sumUpChange_;
674#if INFEAS==1
675    sum += INFEAS_MULTIPLIER*numberTimesUpInfeasible_ * CoinMax(distanceToCutoff / (upCost + 1.0e-12), sumUpCost_);
676#endif
677#elif TYPE2==1
678    abort();
679#if INFEAS==1
680    sum += INFEAS_MULTIPLIER*numberTimesUpInfeasible_ * (distanceToCutoff / (upCost + 1.0e-12));
681#endif
682#endif
683#if MOD_SHADOW>0
684    if (!upShadowPrice_) {
685        if (number > 0.0)
686            upCost *= sum / number;
687        else
688            upCost  *=  upDynamicPseudoCost_;
689    } else if (upShadowPrice_ > 0.0) {
690        upCost *= upShadowPrice_;
691    } else {
692        upCost *= (upDynamicPseudoCost_ - upShadowPrice_);
693    }
694#else
695    if (upShadowPrice_ <= 0.0) {
696        if (number > 0.0)
697            upCost *= sum / number;
698        else
699            upCost  *=  upDynamicPseudoCost_;
700    } else {
701        upCost *= upShadowPrice_;
702    }
703#endif
704    if (downCost >= upCost)
705        preferredWay = 1;
706    else
707        preferredWay = -1;
708    // See if up down choice set
709    if (upDownSeparator_ > 0.0) {
710        preferredWay = (value - below >= upDownSeparator_) ? 1 : -1;
711    }
712#ifdef FUNNY_BRANCHING2
713    if (fabs(value - nearest) > integerTolerance) {
714        double ratio = (100.0 + lastUpDecrease_) / (100.0 + lastDownDecrease_);
715        downCost *= ratio;
716        upCost *= ratio;
717        if ((lastUpDecrease_ % 100) == -1)
718            printf("col %d total %d djtimes %d\n",
719                   columnNumber_, lastDownDecrease_, lastUpDecrease_);
720    }
721#endif
722    if (preferredWay_)
723        preferredWay = preferredWay_;
724    if (info->hotstartSolution_) {
725        double targetValue = info->hotstartSolution_[columnNumber_];
726        if (value > targetValue)
727            preferredWay = -1;
728        else
729            preferredWay = 1;
730    }
731    if (fabs(value - nearest) <= integerTolerance) {
732        if (priority_ != -999)
733            return 0.0;
734        else
735            return 1.0e-13;
736    } else {
737        int stateOfSearch = model_->stateOfSearch() % 10;
738        double returnValue = 0.0;
739        double minValue = CoinMin(downCost, upCost);
740        double maxValue = CoinMax(downCost, upCost);
741#ifdef COIN_DEVELOP
742        char where;
743#endif
744        // was <= 10
745        //if (stateOfSearch<=1||model_->currentNode()->depth()<=-10 /* was ||maxValue>0.2*distanceToCutoff*/) {
746        if (stateOfSearch <= 2) {
747            // no branching solution
748#ifdef COIN_DEVELOP
749            where = 'i';
750#endif
751            returnValue = WEIGHT_BEFORE * minValue + (1.0 - WEIGHT_BEFORE) * maxValue;
752            if (0) {
753                double sum;
754                int number;
755                double downCost2 = CoinMax(value - below, 0.0);
756                sum = sumDownCost_;
757                number = numberTimesDown_;
758                if (number > 0)
759                    downCost2 *= sum / static_cast<double> (number);
760                else
761                    downCost2  *=  downDynamicPseudoCost_;
762                double upCost2 = CoinMax((above - value), 0.0);
763                sum = sumUpCost_;
764                number = numberTimesUp_;
765                if (number > 0)
766                    upCost2 *= sum / static_cast<double> (number);
767                else
768                    upCost2  *=  upDynamicPseudoCost_;
769                double minValue2 = CoinMin(downCost2, upCost2);
770                double maxValue2 = CoinMax(downCost2, upCost2);
771                printf("%d value %g downC %g upC %g minV %g maxV %g downC2 %g upC2 %g minV2 %g maxV2 %g\n",
772                       columnNumber_, value, downCost, upCost, minValue, maxValue,
773                       downCost2, upCost2, minValue2, maxValue2);
774            }
775        } else {
776            // some solution
777#ifdef COIN_DEVELOP
778            where = 'I';
779#endif
780#ifndef WEIGHT_PRODUCT
781            returnValue = WEIGHT_AFTER * minValue + (1.0 - WEIGHT_AFTER) * maxValue;
782#else
783            double minProductWeight = model_->getDblParam(CbcModel::CbcSmallChange);
784            returnValue = CoinMax(minValue, minProductWeight) * CoinMax(maxValue, minProductWeight);
785            //returnValue += minProductWeight*minValue;
786#endif
787        }
788        if (numberTimesUp_ < numberBeforeTrust_ ||
789                numberTimesDown_ < numberBeforeTrust_) {
790            //if (returnValue<1.0e10)
791            //returnValue += 1.0e12;
792            //else
793            returnValue *= 1.0e3;
794            if (!numberTimesUp_ && !numberTimesDown_)
795                returnValue *= 1.0e10;
796        }
797        //if (fabs(value-0.5)<1.0e-5) {
798        //returnValue = 3.0*returnValue + 0.2;
799        //} else if (value>0.9) {
800        //returnValue = 2.0*returnValue + 0.1;
801        //}
802        if (method_ == 1) {
803            // probing
804            // average
805            double up = 1.0e-15;
806            double down = 1.0e-15;
807            if (numberTimesProbingTotal_) {
808                up += numberTimesUpTotalFixed_ / static_cast<double> (numberTimesProbingTotal_);
809                down += numberTimesDownTotalFixed_ / static_cast<double> (numberTimesProbingTotal_);
810            }
811            returnValue = 1 + 10.0 * CoinMin(numberTimesDownLocalFixed_, numberTimesUpLocalFixed_) +
812                          CoinMin(down, up);
813            returnValue *= 1.0e-3;
814        }
815#ifdef COIN_DEVELOP
816        History hist;
817        hist.where_ = where;
818        hist.status_ = ' ';
819        hist.sequence_ = columnNumber_;
820        hist.numberUp_ = numberTimesUp_;
821        hist.numberUpInf_ = numberTimesUpInfeasible_;
822        hist.sumUp_ = sumUpCost_;
823        hist.upEst_ = upCost;
824        hist.numberDown_ = numberTimesDown_;
825        hist.numberDownInf_ = numberTimesDownInfeasible_;
826        hist.sumDown_ = sumDownCost_;
827        hist.downEst_ = downCost;
828        if (stateOfSearch)
829            addRecord(hist);
830#endif
831        return CoinMax(returnValue, 1.0e-15);
832    }
833}
834// Creates a branching object
835CbcBranchingObject *
836CbcSimpleIntegerDynamicPseudoCost::createCbcBranch(OsiSolverInterface * /*solver*/,
837        const OsiBranchingInformation * info, int way)
838{
839    double value = info->solution_[columnNumber_];
840    value = CoinMax(value, info->lower_[columnNumber_]);
841    value = CoinMin(value, info->upper_[columnNumber_]);
842    assert (info->upper_[columnNumber_] > info->lower_[columnNumber_]);
843    if (!info->hotstartSolution_ && priority_ != -999) {
844#ifndef NDEBUG
845#ifndef SWITCH_VARIABLES
846        double nearest = floor(value + 0.5);
847        assert (fabs(value - nearest) > info->integerTolerance_);
848#endif
849#endif
850    } else if (info->hotstartSolution_) {
851        double targetValue = info->hotstartSolution_[columnNumber_];
852        if (way > 0)
853            value = targetValue - 0.1;
854        else
855            value = targetValue + 0.1;
856    } else {
857        if (value <= info->lower_[columnNumber_])
858            value += 0.1;
859        else if (value >= info->upper_[columnNumber_])
860            value -= 0.1;
861    }
862    assert (value >= info->lower_[columnNumber_] &&
863            value <= info->upper_[columnNumber_]);
864    CbcDynamicPseudoCostBranchingObject * newObject =
865        new CbcDynamicPseudoCostBranchingObject(model_, columnNumber_, way,
866                                                value, this);
867    double up =  upDynamicPseudoCost_ * (ceil(value) - value);
868    double down =  downDynamicPseudoCost_ * (value - floor(value));
869    double changeInGuessed = up - down;
870    if (way > 0)
871        changeInGuessed = - changeInGuessed;
872    changeInGuessed = CoinMax(0.0, changeInGuessed);
873    //if (way>0)
874    //changeInGuessed += 1.0e8; // bias to stay up
875    newObject->setChangeInGuessed(changeInGuessed);
876    newObject->setOriginalObject(this);
877    return newObject;
878}
879
880// Return "up" estimate
881double
882CbcSimpleIntegerDynamicPseudoCost::upEstimate() const
883{
884    const double * solution = model_->testSolution();
885    const double * lower = model_->getCbcColLower();
886    const double * upper = model_->getCbcColUpper();
887    double value = solution[columnNumber_];
888    value = CoinMax(value, lower[columnNumber_]);
889    value = CoinMin(value, upper[columnNumber_]);
890    if (upper[columnNumber_] == lower[columnNumber_]) {
891        // fixed
892        return 0.0;
893    }
894    double integerTolerance =
895        model_->getDblParam(CbcModel::CbcIntegerTolerance);
896    double below = floor(value + integerTolerance);
897    double above = below + 1.0;
898    if (above > upper[columnNumber_]) {
899        above = below;
900        below = above - 1;
901    }
902    double upCost = CoinMax((above - value) * upDynamicPseudoCost_, 0.0);
903    return upCost;
904}
905// Return "down" estimate
906double
907CbcSimpleIntegerDynamicPseudoCost::downEstimate() const
908{
909    const double * solution = model_->testSolution();
910    const double * lower = model_->getCbcColLower();
911    const double * upper = model_->getCbcColUpper();
912    double value = solution[columnNumber_];
913    value = CoinMax(value, lower[columnNumber_]);
914    value = CoinMin(value, upper[columnNumber_]);
915    if (upper[columnNumber_] == lower[columnNumber_]) {
916        // fixed
917        return 0.0;
918    }
919    double integerTolerance =
920        model_->getDblParam(CbcModel::CbcIntegerTolerance);
921    double below = floor(value + integerTolerance);
922    double above = below + 1.0;
923    if (above > upper[columnNumber_]) {
924        above = below;
925        below = above - 1;
926    }
927    double downCost = CoinMax((value - below) * downDynamicPseudoCost_, 0.0);
928    return downCost;
929}
930// Set down pseudo cost
931void
932CbcSimpleIntegerDynamicPseudoCost::setDownDynamicPseudoCost(double value)
933{
934#ifdef TRACE_ONE
935    double oldDown = sumDownCost_;
936#endif
937    downDynamicPseudoCost_ = value;
938    sumDownCost_ = CoinMax(sumDownCost_, value * numberTimesDown_);
939#ifdef TRACE_ONE
940    if (columnNumber_ == TRACE_ONE) {
941        double down = downDynamicPseudoCost_ * numberTimesDown_;
942        printf("For %d sumDown %g (%d), inf (%d) - pseudo %g - sumDown was %g -> %g\n",
943               TRACE_ONE, down, numberTimesDown_,
944               numberTimesDownInfeasible_, downDynamicPseudoCost_,
945               oldDown, sumDownCost_);
946    }
947#endif
948}
949// Modify down pseudo cost in a slightly different way
950void
951CbcSimpleIntegerDynamicPseudoCost::updateDownDynamicPseudoCost(double value)
952{
953    sumDownCost_ += value;
954    numberTimesDown_++;
955    downDynamicPseudoCost_ = sumDownCost_ / static_cast<double>(numberTimesDown_);
956}
957// Set up pseudo cost
958void
959CbcSimpleIntegerDynamicPseudoCost::setUpDynamicPseudoCost(double value)
960{
961#ifdef TRACE_ONE
962    double oldUp = sumUpCost_;
963#endif
964    upDynamicPseudoCost_ = value;
965    sumUpCost_ = CoinMax(sumUpCost_, value * numberTimesUp_);
966#ifdef TRACE_ONE
967    if (columnNumber_ == TRACE_ONE) {
968        double up = upDynamicPseudoCost_ * numberTimesUp_;
969        printf("For %d sumUp %g (%d), inf (%d) - pseudo %g - sumUp was %g -> %g\n",
970               TRACE_ONE, up, numberTimesUp_,
971               numberTimesUpInfeasible_, upDynamicPseudoCost_,
972               oldUp, sumUpCost_);
973    }
974#endif
975}
976// Modify up pseudo cost in a slightly different way
977void
978CbcSimpleIntegerDynamicPseudoCost::updateUpDynamicPseudoCost(double value)
979{
980    sumUpCost_ += value;
981    numberTimesUp_++;
982    upDynamicPseudoCost_ = sumUpCost_ / static_cast<double>(numberTimesUp_);
983}
984/* Pass in information on branch just done and create CbcObjectUpdateData instance.
985   If object does not need data then backward pointer will be NULL.
986   Assumes can get information from solver */
987CbcObjectUpdateData
988CbcSimpleIntegerDynamicPseudoCost::createUpdateInformation(const OsiSolverInterface * solver,
989        const CbcNode * node,
990        const CbcBranchingObject * branchingObject)
991{
992    double originalValue = node->objectiveValue();
993    int originalUnsatisfied = node->numberUnsatisfied();
994    double objectiveValue = solver->getObjValue() * solver->getObjSense();
995    int unsatisfied = 0;
996    int i;
997    //might be base model - doesn't matter
998    int numberIntegers = model_->numberIntegers();;
999    const double * solution = solver->getColSolution();
1000    //const double * lower = solver->getColLower();
1001    //const double * upper = solver->getColUpper();
1002    double change = CoinMax(0.0, objectiveValue - originalValue);
1003    int iStatus;
1004    if (solver->isProvenOptimal())
1005        iStatus = 0; // optimal
1006    else if (solver->isIterationLimitReached()
1007             && !solver->isDualObjectiveLimitReached())
1008        iStatus = 2; // unknown
1009    else
1010        iStatus = 1; // infeasible
1011
1012    bool feasible = iStatus != 1;
1013    if (feasible) {
1014        double integerTolerance =
1015            model_->getDblParam(CbcModel::CbcIntegerTolerance);
1016        const int * integerVariable = model_->integerVariable();
1017        for (i = 0; i < numberIntegers; i++) {
1018            int j = integerVariable[i];
1019            double value = solution[j];
1020            double nearest = floor(value + 0.5);
1021            if (fabs(value - nearest) > integerTolerance)
1022                unsatisfied++;
1023#ifdef SWITCH_VARIABLES
1024            const CbcSwitchingBinary * sObject = dynamic_cast<const CbcSwitchingBinary *> (this);
1025            if (sObject) {
1026              int state[3],nBadFixed;
1027              unsatisfied +=
1028                sObject->checkAssociatedBounds(solver,solution,0,
1029                                               state,nBadFixed);
1030            }
1031#endif
1032        }
1033    }
1034    int way = branchingObject->way();
1035    way = - way; // because after branch so moved on
1036    double value = branchingObject->value();
1037    CbcObjectUpdateData newData (this, way,
1038                                 change, iStatus,
1039                                 originalUnsatisfied - unsatisfied, value);
1040    newData.originalObjective_ = originalValue;
1041    // Solvers know about direction
1042    double direction = solver->getObjSense();
1043    solver->getDblParam(OsiDualObjectiveLimit, newData.cutoff_);
1044    newData.cutoff_ *= direction;
1045    return newData;
1046}
1047// Just update using feasible branches and keep count of infeasible
1048#undef INFEAS
1049// Update object by CbcObjectUpdateData
1050void
1051CbcSimpleIntegerDynamicPseudoCost::updateInformation(const CbcObjectUpdateData & data)
1052{
1053    bool feasible = data.status_ != 1;
1054    int way = data.way_;
1055    double value = data.branchingValue_;
1056    double change = data.change_;
1057#ifdef COIN_DEVELOP
1058    History hist;
1059    hist.where_ = 'U'; // need to tell if hot
1060#endif
1061    double movement = 0.0;
1062    if (way < 0) {
1063        // down
1064        movement = value - floor(value);
1065        if (feasible) {
1066#ifdef COIN_DEVELOP
1067            hist.status_ = 'D';
1068#endif
1069            movement = CoinMax(movement, MINIMUM_MOVEMENT);
1070            //printf("(down change %g value down %g ",change,movement);
1071            incrementNumberTimesDown();
1072            addToSumDownChange(1.0e-30 + movement);
1073            addToSumDownDecrease(data.intDecrease_);
1074#if TYPE2==0
1075            addToSumDownCost(change / (1.0e-30 + movement));
1076            setDownDynamicPseudoCost(sumDownCost() / static_cast<double>( numberTimesDown()));
1077#elif TYPE2==1
1078            addToSumDownCost(change);
1079            setDownDynamicPseudoCost(sumDownCost() / sumDownChange());
1080#elif TYPE2==2
1081            addToSumDownCost(change*TYPERATIO + (1.0 - TYPERATIO)*change / (1.0e-30 + movement));
1082            setDownDynamicPseudoCost(sumDownCost()*(TYPERATIO / sumDownChange() + (1.0 - TYPERATIO) / (double) numberTimesDown()));
1083#endif
1084        } else {
1085#ifdef COIN_DEVELOP
1086            hist.status_ = 'd';
1087#endif
1088            //printf("(down infeasible value down %g ",change,movement);
1089            incrementNumberTimesDown();
1090            incrementNumberTimesDownInfeasible();
1091#if INFEAS==2
1092            double distanceToCutoff = 0.0;
1093            double objectiveValue = model->getCurrentMinimizationObjValue();
1094            distanceToCutoff =  model->getCutoff()  - originalValue;
1095            if (distanceToCutoff < 1.0e20)
1096                change = distanceToCutoff * 2.0;
1097            else
1098                change = downDynamicPseudoCost() * movement * 10.0;
1099            change = CoinMax(1.0e-12 * (1.0 + fabs(originalValue)), change);
1100            addToSumDownChange(1.0e-30 + movement);
1101            addToSumDownDecrease(data.intDecrease_);
1102#if TYPE2==0
1103            addToSumDownCost(change / (1.0e-30 + movement));
1104            setDownDynamicPseudoCost(sumDownCost() / (double) numberTimesDown());
1105#elif TYPE2==1
1106            addToSumDownCost(change);
1107            setDownDynamicPseudoCost(sumDownCost() / sumDownChange());
1108#elif TYPE2==2
1109            addToSumDownCost(change*TYPERATIO + (1.0 - TYPERATIO)*change / (1.0e-30 + movement));
1110            setDownDynamicPseudoCost(sumDownCost()*(TYPERATIO / sumDownChange() + (1.0 - TYPERATIO) / (double) numberTimesDown()));
1111#endif
1112#endif
1113        }
1114#if INFEAS==1
1115        double sum = sumDownCost_;
1116        int number = numberTimesDown_;
1117        double originalValue = data.originalObjective_;
1118        assert (originalValue != COIN_DBL_MAX);
1119        double distanceToCutoff =  data.cutoff_  - originalValue;
1120        if (distanceToCutoff > 1.0e20)
1121            distanceToCutoff = 10.0 + fabs(originalValue);
1122        sum += INFEAS_MULTIPLIER*numberTimesDownInfeasible_ * CoinMax(distanceToCutoff, 1.0e-12 * (1.0 + fabs(originalValue)));
1123        setDownDynamicPseudoCost(sum / static_cast<double> (number));
1124#endif
1125    } else {
1126        // up
1127        movement = ceil(value) - value;
1128        if (feasible) {
1129#ifdef COIN_DEVELOP
1130            hist.status_ = 'U';
1131#endif
1132            movement = CoinMax(movement, MINIMUM_MOVEMENT);
1133            //printf("(up change %g value down %g ",change,movement);
1134            incrementNumberTimesUp();
1135            addToSumUpChange(1.0e-30 + movement);
1136            addToSumUpDecrease(data.intDecrease_);
1137#if TYPE2==0
1138            addToSumUpCost(change / (1.0e-30 + movement));
1139            setUpDynamicPseudoCost(sumUpCost() / static_cast<double> (numberTimesUp()));
1140#elif TYPE2==1
1141            addToSumUpCost(change);
1142            setUpDynamicPseudoCost(sumUpCost() / sumUpChange());
1143#elif TYPE2==2
1144            addToSumUpCost(change*TYPERATIO + (1.0 - TYPERATIO)*change / (1.0e-30 + movement));
1145            setUpDynamicPseudoCost(sumUpCost()*(TYPERATIO / sumUpChange() + (1.0 - TYPERATIO) / (double) numberTimesUp()));
1146#endif
1147        } else {
1148#ifdef COIN_DEVELOP
1149            hist.status_ = 'u';
1150#endif
1151            //printf("(up infeasible value down %g ",change,movement);
1152            incrementNumberTimesUp();
1153            incrementNumberTimesUpInfeasible();
1154#if INFEAS==2
1155            double distanceToCutoff = 0.0;
1156            double objectiveValue = model->getCurrentMinimizationObjValue();
1157            distanceToCutoff =  model->getCutoff()  - originalValue;
1158            if (distanceToCutoff < 1.0e20)
1159                change = distanceToCutoff * 2.0;
1160            else
1161                change = upDynamicPseudoCost() * movement * 10.0;
1162            change = CoinMax(1.0e-12 * (1.0 + fabs(originalValue)), change);
1163            addToSumUpChange(1.0e-30 + movement);
1164            addToSumUpDecrease(data.intDecrease_);
1165#if TYPE2==0
1166            addToSumUpCost(change / (1.0e-30 + movement));
1167            setUpDynamicPseudoCost(sumUpCost() / (double) numberTimesUp());
1168#elif TYPE2==1
1169            addToSumUpCost(change);
1170            setUpDynamicPseudoCost(sumUpCost() / sumUpChange());
1171#elif TYPE2==2
1172            addToSumUpCost(change*TYPERATIO + (1.0 - TYPERATIO)*change / (1.0e-30 + movement));
1173            setUpDynamicPseudoCost(sumUpCost()*(TYPERATIO / sumUpChange() + (1.0 - TYPERATIO) / (double) numberTimesUp()));
1174#endif
1175#endif
1176        }
1177#if INFEAS==1
1178        double sum = sumUpCost_;
1179        int number = numberTimesUp_;
1180        double originalValue = data.originalObjective_;
1181        assert (originalValue != COIN_DBL_MAX);
1182        double distanceToCutoff =  data.cutoff_  - originalValue;
1183        if (distanceToCutoff > 1.0e20)
1184            distanceToCutoff = 10.0 + fabs(originalValue);
1185        sum += INFEAS_MULTIPLIER*numberTimesUpInfeasible_ * CoinMax(distanceToCutoff, 1.0e-12 * (1.0 + fabs(originalValue)));
1186        setUpDynamicPseudoCost(sum / static_cast<double> (number));
1187#endif
1188    }
1189    if (data.way_ < 0)
1190        assert (numberTimesDown_ > 0);
1191    else
1192        assert (numberTimesUp_ > 0);
1193    assert (downDynamicPseudoCost_ >= 0.0 && downDynamicPseudoCost_ < 1.0e100);
1194    downDynamicPseudoCost_ = CoinMax(1.0e-10, downDynamicPseudoCost_);
1195    assert (upDynamicPseudoCost_ >= 0.0 && upDynamicPseudoCost_ < 1.0e100);
1196    upDynamicPseudoCost_ = CoinMax(1.0e-10, upDynamicPseudoCost_);
1197#ifdef COIN_DEVELOP
1198    hist.sequence_ = columnNumber_;
1199    hist.numberUp_ = numberTimesUp_;
1200    hist.numberUpInf_ = numberTimesUpInfeasible_;
1201    hist.sumUp_ = sumUpCost_;
1202    hist.upEst_ = change;
1203    hist.numberDown_ = numberTimesDown_;
1204    hist.numberDownInf_ = numberTimesDownInfeasible_;
1205    hist.sumDown_ = sumDownCost_;
1206    hist.downEst_ = movement;
1207    addRecord(hist);
1208#endif
1209    //print(1,0.5);
1210    assert (downDynamicPseudoCost_ > 1.0e-40 && upDynamicPseudoCost_ > 1.0e-40);
1211#if MOD_SHADOW>1
1212    if (upShadowPrice_ > 0.0 && numberTimesDown_ >= numberBeforeTrust_
1213            && numberTimesUp_ >= numberBeforeTrust_) {
1214        // Set negative
1215        upShadowPrice_ = -upShadowPrice_;
1216        assert (downShadowPrice_ > 0.0);
1217        downShadowPrice_ = - downShadowPrice_;
1218    }
1219#endif
1220}
1221// Updates stuff like pseudocosts after mini branch and bound
1222void
1223CbcSimpleIntegerDynamicPseudoCost::updateAfterMini(int numberDown, int numberDownInfeasible,
1224        double sumDown, int numberUp,
1225        int numberUpInfeasible, double sumUp)
1226{
1227    numberTimesDown_ = numberDown;
1228    numberTimesDownInfeasible_ = numberDownInfeasible;
1229    sumDownCost_ = sumDown;
1230    numberTimesUp_ = numberUp;
1231    numberTimesUpInfeasible_ = numberUpInfeasible;
1232    sumUpCost_ = sumUp;
1233    if (numberTimesDown_ > 0) {
1234        setDownDynamicPseudoCost(sumDownCost_ / static_cast<double> (numberTimesDown_));
1235        assert (downDynamicPseudoCost_ > 0.0 && downDynamicPseudoCost_ < 1.0e50);
1236    }
1237    if (numberTimesUp_ > 0) {
1238        setUpDynamicPseudoCost(sumUpCost_ / static_cast<double> (numberTimesUp_));
1239        assert (upDynamicPseudoCost_ > 0.0 && upDynamicPseudoCost_ < 1.0e50);
1240    }
1241    assert (downDynamicPseudoCost_ > 1.0e-40 && upDynamicPseudoCost_ > 1.0e-40);
1242}
1243// Pass in probing information
1244void
1245CbcSimpleIntegerDynamicPseudoCost::setProbingInformation(int fixedDown, int fixedUp)
1246{
1247    numberTimesProbingTotal_++;
1248    numberTimesDownLocalFixed_ = fixedDown;
1249    numberTimesDownTotalFixed_ += fixedDown;
1250    numberTimesUpLocalFixed_ = fixedUp;
1251    numberTimesUpTotalFixed_ += fixedUp;
1252}
1253// Print
1254void
1255CbcSimpleIntegerDynamicPseudoCost::print(int type, double value) const
1256{
1257    if (!type) {
1258        double meanDown = 0.0;
1259        double devDown = 0.0;
1260        if (numberTimesDown_) {
1261            meanDown = sumDownCost_ / static_cast<double> (numberTimesDown_);
1262            devDown = meanDown * meanDown - 2.0 * meanDown * sumDownCost_;
1263            if (devDown >= 0.0)
1264                devDown = sqrt(devDown);
1265        }
1266        double meanUp = 0.0;
1267        double devUp = 0.0;
1268        if (numberTimesUp_) {
1269            meanUp = sumUpCost_ / static_cast<double> (numberTimesUp_);
1270            devUp = meanUp * meanUp - 2.0 * meanUp * sumUpCost_;
1271            if (devUp >= 0.0)
1272                devUp = sqrt(devUp);
1273        }
1274        printf("%d down %d times (%d inf) mean %g (dev %g) up %d times (%d inf) mean %g (dev %g)\n",
1275               columnNumber_,
1276               numberTimesDown_, numberTimesDownInfeasible_, meanDown, devDown,
1277               numberTimesUp_, numberTimesUpInfeasible_, meanUp, devUp);
1278    } else {
1279        const double * upper = model_->getCbcColUpper();
1280        double integerTolerance =
1281            model_->getDblParam(CbcModel::CbcIntegerTolerance);
1282        double below = floor(value + integerTolerance);
1283        double above = below + 1.0;
1284        if (above > upper[columnNumber_]) {
1285            above = below;
1286            below = above - 1;
1287        }
1288        double objectiveValue = model_->getCurrentMinimizationObjValue();
1289        double distanceToCutoff =  model_->getCutoff() - objectiveValue;
1290        if (distanceToCutoff < 1.0e20)
1291            distanceToCutoff *= 10.0;
1292        else
1293            distanceToCutoff = 1.0e2 + fabs(objectiveValue);
1294        distanceToCutoff = CoinMax(distanceToCutoff, 1.0e-12 * (1.0 + fabs(objectiveValue)));
1295        double sum;
1296        int number;
1297        double downCost = CoinMax(value - below, 0.0);
1298        double downCost0 = downCost * downDynamicPseudoCost_;
1299        sum = sumDownCost();
1300        number = numberTimesDown();
1301        sum += INFEAS_MULTIPLIER*numberTimesDownInfeasible() * (distanceToCutoff / (downCost + 1.0e-12));
1302        if (number > 0)
1303            downCost *= sum / static_cast<double> (number);
1304        else
1305            downCost  *=  downDynamicPseudoCost_;
1306        double upCost = CoinMax((above - value), 0.0);
1307        double upCost0 = upCost * upDynamicPseudoCost_;
1308        sum = sumUpCost();
1309        number = numberTimesUp();
1310        sum += INFEAS_MULTIPLIER*numberTimesUpInfeasible() * (distanceToCutoff / (upCost + 1.0e-12));
1311        if (number > 0)
1312            upCost *= sum / static_cast<double> (number);
1313        else
1314            upCost  *=  upDynamicPseudoCost_;
1315        printf("%d down %d times %g (est %g)  up %d times %g (est %g)\n",
1316               columnNumber_,
1317               numberTimesDown_, downCost, downCost0,
1318               numberTimesUp_, upCost, upCost0);
1319    }
1320}
1321
1322//##############################################################################
1323
1324// Default Constructor
1325CbcIntegerPseudoCostBranchingObject::CbcIntegerPseudoCostBranchingObject()
1326        : CbcIntegerBranchingObject()
1327{
1328    changeInGuessed_ = 1.0e-5;
1329}
1330
1331// Useful constructor
1332CbcIntegerPseudoCostBranchingObject::CbcIntegerPseudoCostBranchingObject (CbcModel * model,
1333        int variable, int way , double value)
1334        : CbcIntegerBranchingObject(model, variable, way, value)
1335{
1336}
1337// Useful constructor for fixing
1338CbcIntegerPseudoCostBranchingObject::CbcIntegerPseudoCostBranchingObject (CbcModel * model,
1339        int variable, int way,
1340        double lowerValue,
1341        double /*upperValue*/)
1342        : CbcIntegerBranchingObject(model, variable, way, lowerValue)
1343{
1344    changeInGuessed_ = 1.0e100;
1345}
1346
1347
1348// Copy constructor
1349CbcIntegerPseudoCostBranchingObject::CbcIntegerPseudoCostBranchingObject (
1350    const CbcIntegerPseudoCostBranchingObject & rhs)
1351        : CbcIntegerBranchingObject(rhs)
1352{
1353    changeInGuessed_ = rhs.changeInGuessed_;
1354}
1355
1356// Assignment operator
1357CbcIntegerPseudoCostBranchingObject &
1358CbcIntegerPseudoCostBranchingObject::operator=( const CbcIntegerPseudoCostBranchingObject & rhs)
1359{
1360    if (this != &rhs) {
1361        CbcIntegerBranchingObject::operator=(rhs);
1362        changeInGuessed_ = rhs.changeInGuessed_;
1363    }
1364    return *this;
1365}
1366CbcBranchingObject *
1367CbcIntegerPseudoCostBranchingObject::clone() const
1368{
1369    return (new CbcIntegerPseudoCostBranchingObject(*this));
1370}
1371
1372
1373// Destructor
1374CbcIntegerPseudoCostBranchingObject::~CbcIntegerPseudoCostBranchingObject ()
1375{
1376}
1377
1378/*
1379  Perform a branch by adjusting the bounds of the specified variable. Note
1380  that each arm of the branch advances the object to the next arm by
1381  advancing the value of way_.
1382
1383  Providing new values for the variable's lower and upper bounds for each
1384  branching direction gives a little bit of additional flexibility and will
1385  be easily extensible to multi-way branching.
1386  Returns change in guessed objective on next branch
1387*/
1388double
1389CbcIntegerPseudoCostBranchingObject::branch()
1390{
1391    CbcIntegerBranchingObject::branch();
1392    return changeInGuessed_;
1393}
1394
1395/** Compare the \c this with \c brObj. \c this and \c brObj must be os the
1396    same type and must have the same original object, but they may have
1397    different feasible regions.
1398    Return the appropriate CbcRangeCompare value (first argument being the
1399    sub/superset if that's the case). In case of overlap (and if \c
1400    replaceIfOverlap is true) replace the current branching object with one
1401    whose feasible region is the overlap.
1402*/
1403CbcRangeCompare
1404CbcIntegerPseudoCostBranchingObject::compareBranchingObject
1405(const CbcBranchingObject* brObj, const bool replaceIfOverlap)
1406{
1407    const CbcIntegerPseudoCostBranchingObject* br =
1408        dynamic_cast<const CbcIntegerPseudoCostBranchingObject*>(brObj);
1409    assert(br);
1410    double* thisBd = way_ < 0 ? down_ : up_;
1411    const double* otherBd = br->way_ < 0 ? br->down_ : br->up_;
1412    return CbcCompareRanges(thisBd, otherBd, replaceIfOverlap);
1413}
1414#ifdef SWITCH_VARIABLES
1415/** Default Constructor
1416
1417  Equivalent to an unspecified binary variable.
1418*/
1419CbcSwitchingBinary::CbcSwitchingBinary ()
1420        : CbcSimpleIntegerDynamicPseudoCost(),
1421          zeroLowerBound_(NULL),
1422          oneLowerBound_(NULL),
1423          zeroUpperBound_(NULL),
1424          oneUpperBound_(NULL),
1425          otherVariable_(NULL),
1426          numberOther_(0),
1427          type_(0)
1428{
1429}
1430
1431/** Useful constructor
1432*/
1433CbcSwitchingBinary::CbcSwitchingBinary (CbcSimpleIntegerDynamicPseudoCost * oldObject,
1434                                        int nOdd,const int * other, const int * otherRow)
1435  : CbcSimpleIntegerDynamicPseudoCost(*oldObject),
1436    zeroLowerBound_(NULL),
1437    oneLowerBound_(NULL),
1438    zeroUpperBound_(NULL),
1439    oneUpperBound_(NULL),
1440    otherVariable_(NULL),
1441    numberOther_(0),
1442    type_(0)
1443{
1444  if (nOdd)
1445    type_=2;
1446  const CoinPackedMatrix * rowCopy = model_->solver()->getMatrixByRow();
1447  const int * column = rowCopy->getIndices();
1448  //const int * rowLength = rowCopy->getVectorLengths();
1449  const CoinBigIndex * rowStart = rowCopy->getVectorStarts();
1450  //const double * rowLower = model_->solver()->getRowLower();
1451  const double * rowUpper = model_->solver()->getRowUpper();
1452  const double * columnLower = model_->solver()->getColLower();
1453  const double * columnUpper = model_->solver()->getColUpper();
1454  const double * element = rowCopy->getElements();
1455  int last = other[0];
1456  int nPair=0;
1457  int nInGroup=1;
1458  for (int i=1;i<=nOdd;i++) {
1459    if (other[i]==last) {
1460      nInGroup++;
1461    } else {
1462      if (nInGroup>2 && model_->logLevel()>2)
1463        printf("%d in group for column %d - some redundancy\n",
1464               nInGroup,columnNumber_);
1465      nPair++;
1466      last=other[i];
1467      nInGroup=1;
1468    }
1469  }
1470  zeroLowerBound_ = new double [4*nPair];
1471  oneLowerBound_ = zeroLowerBound_+nPair;
1472  zeroUpperBound_ = oneLowerBound_+nPair;
1473  oneUpperBound_ = zeroUpperBound_+nPair;
1474  otherVariable_ = new int [nPair];
1475  numberOther_ = nPair;
1476  if (nPair>1&&model_->logLevel()>2)
1477    printf("%d pairs for column %d\n",
1478               nPair,columnNumber_);
1479  // Now fill
1480  last = other[0];
1481  nPair=0;
1482  int rows[20];
1483  rows[0]= otherRow[0];
1484  nInGroup=1;
1485  for (int i=1;i<=nOdd;i++) {
1486    if (other[i]==last) {
1487      rows[nInGroup++]=otherRow[i];
1488    } else {
1489      double newLowerZero=0.0;
1490      double newUpperZero=COIN_DBL_MAX;
1491      double newLowerOne=0.0;
1492      double newUpperOne=COIN_DBL_MAX;
1493      int cColumn=-1;
1494      for (int j=0;j<nInGroup;j++) {
1495        int iRow = rows[j];
1496        CoinBigIndex k = rowStart[iRow];
1497        double bValue, cValue;
1498        if (column[k]==columnNumber_) {
1499          bValue=element[k];
1500          cValue=element[k+1];
1501          cColumn=column[k+1];
1502        } else {
1503          bValue=element[k+1];
1504          cValue=element[k];
1505          cColumn=column[k];
1506        }
1507        if (rowUpper[iRow]) {
1508          // G row - convert to L
1509          bValue=-bValue;
1510          cValue=-cValue;
1511        }
1512        if (bValue>0.0) {
1513          // binary*abs(bValue) <= continuous*abs(cValue);
1514          newLowerOne = -bValue/cValue;
1515        } else {
1516          // binary*abs(bValue) >= continuous*abs(cValue);
1517          newUpperOne = -bValue/cValue;
1518          newUpperZero = 0.0;
1519        }
1520      }
1521      zeroLowerBound_[nPair]=newLowerZero;
1522      oneLowerBound_[nPair]=newLowerOne;
1523      zeroUpperBound_[nPair]=newUpperZero;
1524      oneUpperBound_[nPair]=newUpperOne;
1525      // make current bounds tight
1526      double newLower = CoinMin(newLowerZero,newLowerOne);
1527      if (newLower>columnLower[cColumn])
1528        model_->solver()->setColLower(cColumn,newLower);
1529      double newUpper = CoinMax(newUpperZero,newUpperOne);
1530      if (newUpper<columnUpper[cColumn])
1531        model_->solver()->setColUpper(cColumn,newUpper);
1532      otherVariable_[nPair++]=cColumn;
1533      last=other[i];
1534      rows[0] = otherRow[i];
1535      nInGroup=1;
1536    }
1537  }
1538}
1539// Copy constructor
1540CbcSwitchingBinary::CbcSwitchingBinary ( const CbcSwitchingBinary & rhs)
1541        : CbcSimpleIntegerDynamicPseudoCost(rhs),
1542          numberOther_(rhs.numberOther_),
1543          type_(rhs.type_)
1544{
1545  zeroLowerBound_ = CoinCopyOfArray(rhs.zeroLowerBound_,4*numberOther_);
1546  oneLowerBound_ = zeroLowerBound_+numberOther_;
1547  zeroUpperBound_ = oneLowerBound_+numberOther_;
1548  oneUpperBound_ = zeroUpperBound_+numberOther_;
1549  otherVariable_ = CoinCopyOfArray(rhs.otherVariable_,numberOther_);
1550}
1551
1552// Clone
1553CbcObject *
1554CbcSwitchingBinary::clone() const
1555{
1556    return new CbcSwitchingBinary(*this);
1557}
1558
1559// Assignment operator
1560CbcSwitchingBinary &
1561CbcSwitchingBinary::operator=( const CbcSwitchingBinary & rhs)
1562{
1563    if (this != &rhs) {
1564        CbcSimpleIntegerDynamicPseudoCost::operator=(rhs);
1565        numberOther_=rhs.numberOther_;
1566        type_ = rhs.type_;
1567        delete [] zeroLowerBound_;
1568        delete [] otherVariable_;
1569        zeroLowerBound_ = CoinCopyOfArray(rhs.zeroLowerBound_,4*numberOther_);
1570        oneLowerBound_ = zeroLowerBound_+numberOther_;
1571        zeroUpperBound_ = oneLowerBound_+numberOther_;
1572        oneUpperBound_ = zeroUpperBound_+numberOther_;
1573        otherVariable_ = CoinCopyOfArray(rhs.otherVariable_,numberOther_);
1574    }
1575    return *this;
1576}
1577
1578// Destructor
1579CbcSwitchingBinary::~CbcSwitchingBinary ()
1580{
1581  delete [] zeroLowerBound_;
1582  delete [] otherVariable_;
1583}
1584// Add in zero switches
1585void 
1586CbcSwitchingBinary::addZeroSwitches(int nAdd,const int * columns)
1587{
1588  type_ |= 1;
1589  int nNew = numberOther_+nAdd;
1590  double * bounds = new double[4*nNew];
1591  int * other = new int [nNew];
1592  memcpy(other,otherVariable_,numberOther_*sizeof(int));
1593  delete [] otherVariable_;
1594  otherVariable_=other;
1595  memcpy(bounds,zeroLowerBound_,numberOther_*sizeof(double));
1596  memcpy(bounds+nNew,oneLowerBound_,numberOther_*sizeof(double));
1597  memcpy(bounds+2*nNew,zeroUpperBound_,numberOther_*sizeof(double));
1598  memcpy(bounds+3*nNew,oneUpperBound_,numberOther_*sizeof(double));
1599  delete [] zeroLowerBound_;
1600  zeroLowerBound_ = bounds;
1601  oneLowerBound_ = zeroLowerBound_+nNew;
1602  zeroUpperBound_ = oneLowerBound_+nNew;
1603  oneUpperBound_ = zeroUpperBound_+nNew;
1604  for (int i=0;i<nAdd;i++) {
1605    zeroLowerBound_[numberOther_]=0.0;
1606    oneLowerBound_[numberOther_]=0.0;
1607    zeroUpperBound_[numberOther_]=0.0;
1608    oneUpperBound_[numberOther_]=COIN_DBL_MAX;
1609    otherVariable_[numberOther_++]=columns[i];
1610  }
1611}
1612// Same - returns true if contents match(ish)
1613bool
1614CbcSwitchingBinary::same(const CbcSwitchingBinary * otherObject) const
1615{
1616  bool okay = CbcSimpleIntegerDynamicPseudoCost::same(otherObject);
1617    return okay;
1618}
1619double
1620CbcSwitchingBinary::infeasibility(const OsiBranchingInformation * info,
1621        int &preferredWay) const
1622{
1623    assert (downDynamicPseudoCost_ > 1.0e-40 && upDynamicPseudoCost_ > 1.0e-40);
1624    double * solution = const_cast<double *>(model_->testSolution());
1625    const double * lower = model_->getCbcColLower();
1626    const double * upper = model_->getCbcColUpper();
1627    double saveValue = solution[columnNumber_];
1628    if (!lower[columnNumber_]&&upper[columnNumber_]==1.0) {
1629      double integerTolerance =
1630        model_->getDblParam(CbcModel::CbcIntegerTolerance);
1631      if (saveValue<integerTolerance) {
1632        // check others OK
1633        bool allGood=true;
1634        double tolerance;
1635        model_->solver()->getDblParam(OsiPrimalTolerance, tolerance) ;
1636        for (int i=0;i<numberOther_;i++) {
1637          int otherColumn = otherVariable_[i];
1638          double value = solution[otherColumn];
1639          if (value<zeroLowerBound_[i]-tolerance||
1640              value>zeroUpperBound_[i]+tolerance)
1641            allGood=false;
1642        }
1643        if (!allGood)
1644          solution[columnNumber_]=2.0*integerTolerance;
1645      } else if (saveValue>1.0-integerTolerance) {
1646        // check others OK
1647        bool allGood=true;
1648        double tolerance;
1649        model_->solver()->getDblParam(OsiPrimalTolerance, tolerance) ;
1650        for (int i=0;i<numberOther_;i++) {
1651          int otherColumn = otherVariable_[i];
1652          double value = solution[otherColumn];
1653          if (value<oneLowerBound_[i]-tolerance||
1654              value>oneUpperBound_[i]+tolerance)
1655            allGood=false;
1656        }
1657        if (!allGood)
1658          solution[columnNumber_]=1.0-2.0*integerTolerance;
1659      }
1660    }
1661    double inf = CbcSimpleIntegerDynamicPseudoCost::infeasibility(info,preferredWay);
1662    solution[columnNumber_]=saveValue;
1663    return inf;
1664}
1665// Set associated bounds
1666int 
1667CbcSwitchingBinary::setAssociatedBounds(OsiSolverInterface * solver,
1668                                        int cleanBasis) const
1669{
1670  if (!solver)
1671    solver = model_->solver();
1672#ifdef COIN_HAS_CLP
1673  OsiClpSolverInterface * clpSolver
1674    = dynamic_cast<OsiClpSolverInterface *> (solver);
1675  if (cleanBasis!=1)
1676    clpSolver=NULL;
1677#endif
1678  const double * columnLower = solver->getColLower();
1679  const double * columnUpper = solver->getColUpper();
1680  int nChanged=0;
1681  if (!columnUpper[columnNumber_]) {
1682#ifdef COIN_HAS_CLP
1683    if (clpSolver)
1684      clpSolver->setColumnStatus(columnNumber_,ClpSimplex::isFixed);
1685#endif
1686    for (int i=0;i<numberOther_;i++) {
1687      int otherColumn = otherVariable_[i];
1688      if (zeroLowerBound_[i]>columnLower[otherColumn]) {
1689        solver->setColLower(otherColumn,zeroLowerBound_[i]);
1690        nChanged++;
1691      }
1692      if (zeroUpperBound_[i]<columnUpper[otherColumn]) {
1693        solver->setColUpper(otherColumn,zeroUpperBound_[i]);
1694#ifdef COIN_DEVELOP
1695        const double * solution = solver->getColSolution();
1696        double value = solution[otherColumn];
1697        if (value - zeroUpperBound_[i] > 1.0e-5 && model_->logLevel()>1)
1698          printf("value for continuous %d %g - above %g - switch %d is %.12g (ub 0)\n", 
1699                 otherColumn, value, zeroUpperBound_[i],columnNumber_,solution[columnNumber_]);
1700#endif
1701        nChanged++;
1702      }
1703    }
1704  } else if (columnLower[columnNumber_]==1.0) {
1705#ifdef COIN_HAS_CLP
1706    if (clpSolver)
1707      clpSolver->setColumnStatus(columnNumber_,ClpSimplex::isFixed);
1708#endif
1709    for (int i=0;i<numberOther_;i++) {
1710      int otherColumn = otherVariable_[i];
1711      if (oneLowerBound_[i]>columnLower[otherColumn]) {
1712        solver->setColLower(otherColumn,oneLowerBound_[i]);
1713        nChanged++;
1714      }
1715      if (oneUpperBound_[i]<columnUpper[otherColumn]) {
1716        solver->setColUpper(otherColumn,oneUpperBound_[i]);
1717        nChanged++;
1718      }
1719    }
1720  } else if (cleanBasis>=2) {
1721    // if all OK then can fix
1722    int state[3];
1723    int nBadFixed;
1724    const double * solution = solver->getColSolution();
1725    if (!checkAssociatedBounds(solver,solution,
1726                               0,state,nBadFixed)) {
1727      const double *reducedCost = solver->getReducedCost() ;
1728      double good=true;
1729      double integerTolerance =
1730        model_->getDblParam(CbcModel::CbcIntegerTolerance);
1731      if (solution[columnNumber_]<integerTolerance) {
1732        if (cleanBasis==2||reducedCost[columnNumber_]>1.0e-6)
1733          solver->setColUpper(columnNumber_,0.0);
1734        else
1735          good=false;
1736      } else if (solution[columnNumber_]>1.0-integerTolerance) {
1737        if (cleanBasis==2||reducedCost[columnNumber_]<-1.0e-6)
1738          solver->setColLower(columnNumber_,1.0);
1739        else
1740          good=false;
1741      }
1742      if (good)
1743        nChanged=setAssociatedBounds(solver,0);
1744    }
1745  } else {
1746    // see if any continuous bounds force binary
1747    for (int i=0;i<numberOther_;i++) {
1748      int otherColumn = otherVariable_[i];
1749      if (columnLower[otherColumn]>zeroUpperBound_[i]) {
1750        // can't be zero
1751        solver->setColLower(columnNumber_,1.0);
1752        nChanged++;
1753      } else if (columnLower[otherColumn]>oneUpperBound_[i]) {
1754        // can't be one
1755        solver->setColUpper(columnNumber_,0.0);
1756        nChanged++;
1757      }
1758      if (columnUpper[otherColumn]<zeroLowerBound_[i]) {
1759        // can't be zero
1760        solver->setColLower(columnNumber_,1.0);
1761        nChanged++;
1762      } else if (columnUpper[otherColumn]<oneLowerBound_[i]) {
1763        // can't be one
1764        solver->setColUpper(columnNumber_,0.0);
1765        nChanged++;
1766      }
1767    }
1768  }
1769  return nChanged;
1770}
1771// Check associated bounds
1772int 
1773CbcSwitchingBinary::checkAssociatedBounds(const OsiSolverInterface * solver,
1774                                          const double * solution, int printLevel, 
1775                                          int state[3], int & nBadFixed) const
1776{
1777  state[0] = 0;
1778  int nBad=0;
1779  if (!solver)
1780    solver = model_->solver();
1781  double tolerance;
1782  solver->getDblParam(OsiPrimalTolerance, tolerance) ;
1783  const double * columnLower = solver->getColLower();
1784  const double * columnUpper = solver->getColUpper();
1785  double integerTolerance =
1786    model_->getDblParam(CbcModel::CbcIntegerTolerance);
1787  bool printIt = printLevel>2 && model_->logLevel()>1;
1788  if (solution[columnNumber_]<integerTolerance) {
1789    state[0] = -1;
1790    for (int i=0;i<numberOther_;i++) {
1791      int otherColumn = otherVariable_[i];
1792      if (zeroLowerBound_[i]>solution[otherColumn]+tolerance*5.0) {
1793        nBad++;
1794        if (columnUpper[columnNumber_]==0.0) {
1795          nBadFixed++;
1796          //printIt=true;
1797        }
1798        if (printIt)
1799          printf("switch %d at zero, other %d at %.12g below bound of %.12g\n",
1800                 columnNumber_,otherColumn,solution[otherColumn],zeroLowerBound_[i]);
1801      }
1802      if (zeroUpperBound_[i]<solution[otherColumn]-tolerance*5.0) {
1803        nBad++;
1804        if (columnUpper[columnNumber_]==0.0) {
1805          nBadFixed++;
1806          //printIt=true;
1807        }
1808        if (printIt)
1809          printf("switch %d at zero, other %d at %.12g above bound of %.12g\n",
1810                 columnNumber_,otherColumn,solution[otherColumn],zeroUpperBound_[i]);
1811      }
1812    }
1813  } else if (solution[columnNumber_]>1.0-integerTolerance) {
1814    state[0] = 1;
1815    for (int i=0;i<numberOther_;i++) {
1816      int otherColumn = otherVariable_[i];
1817      if (oneLowerBound_[i]>solution[otherColumn]+tolerance*5.0) {
1818        nBad++;
1819        if (columnLower[columnNumber_]==1.0) {
1820          nBadFixed++;
1821          //printIt=true;
1822        }
1823        if (printIt)
1824          printf("switch %d at one, other %d at %.12g below bound of %.12g\n",
1825                 columnNumber_,otherColumn,solution[otherColumn],oneLowerBound_[i]);
1826      }
1827      if (oneUpperBound_[i]<solution[otherColumn]-tolerance*5.0) {
1828        nBad++;
1829        if (columnLower[columnNumber_]==1.0) {
1830          nBadFixed++;
1831          //printIt=true;
1832        }
1833        if (printIt)
1834          printf("switch %d at one, other %d at %.12g above bound of %.12g\n",
1835                 columnNumber_,otherColumn,solution[otherColumn],oneUpperBound_[i]);
1836      }
1837    }
1838  } else {
1839    // in between - compute tight variables
1840    state[1]=0;
1841    state[2]=0;
1842    // for now just compute ones away from bounds
1843    for (int i=0;i<numberOther_;i++) {
1844      int otherColumn = otherVariable_[i];
1845      double otherValue = solution[otherColumn];
1846      if (otherValue>columnLower[otherColumn]+tolerance&&
1847          otherValue<columnUpper[otherColumn]-tolerance)
1848        state[1]++;
1849    }
1850  }
1851  return nBad;
1852}
1853#endif
Note: See TracBrowser for help on using the repository browser.