source: trunk/Cbc/src/CbcSimpleIntegerDynamicPseudoCost.cpp

Last change on this file was 2467, checked in by unxusr, 11 months ago

spaces after angles

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