1 // $Id$
2 // Copyright (C) 2002, International Business Machines
3 // Corporation and others.  All Rights Reserved.
4 // This code is licensed under the terms of the Eclipse Public License (EPL).
5 
6 // Edwin 11/17/2009 - carved out of CbcBranchDynamic
7 
8 #if defined(_MSC_VER)
9 // Turn off compiler warning about long names
10 #pragma warning(disable : 4786)
11 #endif
12 #include <cassert>
13 #include <cstdlib>
14 #include <cmath>
15 #include <cfloat>
16 //#define CBC_DEBUG
17 //#define TRACE_ONE 19
18 #include "OsiSolverInterface.hpp"
19 #include "OsiSolverBranch.hpp"
20 #include "CbcModel.hpp"
21 #include "CbcMessage.hpp"
22 #include "CbcBranchDynamic.hpp"
23 #include "CoinSort.hpp"
24 #include "CoinError.hpp"
25 #include "CbcSimpleIntegerDynamicPseudoCost.hpp"
26 #ifdef COIN_HAS_CLP
27 #include "OsiClpSolverInterface.hpp"
28 #endif
29 #ifdef COIN_DEVELOP
30 typedef 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;
43 History *history = NULL;
44 int numberHistory = 0;
45 int maxHistory = 0;
46 bool getHistoryStatistics_ = true;
increaseHistory()47 static 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 }
addRecord(History newOne)57 static 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 */
CbcSimpleIntegerDynamicPseudoCost()118 CbcSimpleIntegerDynamicPseudoCost::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 */
CbcSimpleIntegerDynamicPseudoCost(CbcModel * model,int iColumn,double breakEven)153 CbcSimpleIntegerDynamicPseudoCost::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 */
CbcSimpleIntegerDynamicPseudoCost(CbcModel * model,int iColumn,double downDynamicPseudoCost,double upDynamicPseudoCost)216 CbcSimpleIntegerDynamicPseudoCost::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 */
CbcSimpleIntegerDynamicPseudoCost(CbcModel * model,int,int iColumn,double downDynamicPseudoCost,double upDynamicPseudoCost)278 CbcSimpleIntegerDynamicPseudoCost::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
CbcSimpleIntegerDynamicPseudoCost(const CbcSimpleIntegerDynamicPseudoCost & rhs)287 CbcSimpleIntegerDynamicPseudoCost::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
320 CbcObject *
clone() const321 CbcSimpleIntegerDynamicPseudoCost::clone() const
322 {
323   return new CbcSimpleIntegerDynamicPseudoCost(*this);
324 }
325 
326 // Assignment operator
327 CbcSimpleIntegerDynamicPseudoCost &
operator =(const CbcSimpleIntegerDynamicPseudoCost & rhs)328 CbcSimpleIntegerDynamicPseudoCost::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
~CbcSimpleIntegerDynamicPseudoCost()363 CbcSimpleIntegerDynamicPseudoCost::~CbcSimpleIntegerDynamicPseudoCost()
364 {
365 }
366 // Copy some information i.e. just variable stuff
copySome(const CbcSimpleIntegerDynamicPseudoCost * otherObject)367 void 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
updateBefore(const OsiObject * rhs)394 void 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
updateAfter(const OsiObject * rhs,const OsiObject * baseObjectX)405 void 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)
same(const CbcSimpleIntegerDynamicPseudoCost * otherObject) const463 bool 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 
516 This returns NULL if branch not represented by bound changes
517 */
518 OsiSolverBranch *
solverBranch() const519 CbcSimpleIntegerDynamicPseudoCost::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
539 double
infeasibility(const OsiBranchingInformation * info,int & preferredWay) const540 CbcSimpleIntegerDynamicPseudoCost::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
821 CbcBranchingObject *
createCbcBranch(OsiSolverInterface *,const OsiBranchingInformation * info,int way)822 CbcSimpleIntegerDynamicPseudoCost::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
865 double
upEstimate() const866 CbcSimpleIntegerDynamicPseudoCost::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
889 double
downEstimate() const890 CbcSimpleIntegerDynamicPseudoCost::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
setDownDynamicPseudoCost(double value)913 void 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
updateDownDynamicPseudoCost(double value)931 void CbcSimpleIntegerDynamicPseudoCost::updateDownDynamicPseudoCost(double value)
932 {
933   sumDownCost_ += value;
934   numberTimesDown_++;
935   downDynamicPseudoCost_ = sumDownCost_ / static_cast< double >(numberTimesDown_);
936 }
937 // Set up pseudo cost
setUpDynamicPseudoCost(double value)938 void 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
updateUpDynamicPseudoCost(double value)956 void 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 */
965 CbcObjectUpdateData
createUpdateInformation(const OsiSolverInterface * solver,const CbcNode * node,const CbcBranchingObject * branchingObject)966 CbcSimpleIntegerDynamicPseudoCost::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
updateInformation(const CbcObjectUpdateData & data)1027 void 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
updateAfterMini(int numberDown,int numberDownInfeasible,double sumDown,int numberUp,int numberUpInfeasible,double sumUp)1198 void 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
setProbingInformation(int fixedDown,int fixedUp)1219 void CbcSimpleIntegerDynamicPseudoCost::setProbingInformation(int fixedDown, int fixedUp)
1220 {
1221   numberTimesProbingTotal_++;
1222   numberTimesDownLocalFixed_ = fixedDown;
1223   numberTimesDownTotalFixed_ += fixedDown;
1224   numberTimesUpLocalFixed_ = fixedUp;
1225   numberTimesUpTotalFixed_ += fixedUp;
1226 }
1227 // Print
print(int type,double value) const1228 void 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
CbcIntegerPseudoCostBranchingObject()1297 CbcIntegerPseudoCostBranchingObject::CbcIntegerPseudoCostBranchingObject()
1298   : CbcIntegerBranchingObject()
1299 {
1300   changeInGuessed_ = 1.0e-5;
1301 }
1302 
1303 // Useful constructor
CbcIntegerPseudoCostBranchingObject(CbcModel * model,int variable,int way,double value)1304 CbcIntegerPseudoCostBranchingObject::CbcIntegerPseudoCostBranchingObject(CbcModel *model,
1305   int variable, int way, double value)
1306   : CbcIntegerBranchingObject(model, variable, way, value)
1307 {
1308 }
1309 // Useful constructor for fixing
CbcIntegerPseudoCostBranchingObject(CbcModel * model,int variable,int way,double lowerValue,double)1310 CbcIntegerPseudoCostBranchingObject::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
CbcIntegerPseudoCostBranchingObject(const CbcIntegerPseudoCostBranchingObject & rhs)1320 CbcIntegerPseudoCostBranchingObject::CbcIntegerPseudoCostBranchingObject(
1321   const CbcIntegerPseudoCostBranchingObject &rhs)
1322   : CbcIntegerBranchingObject(rhs)
1323 {
1324   changeInGuessed_ = rhs.changeInGuessed_;
1325 }
1326 
1327 // Assignment operator
1328 CbcIntegerPseudoCostBranchingObject &
operator =(const CbcIntegerPseudoCostBranchingObject & rhs)1329 CbcIntegerPseudoCostBranchingObject::operator=(const CbcIntegerPseudoCostBranchingObject &rhs)
1330 {
1331   if (this != &rhs) {
1332     CbcIntegerBranchingObject::operator=(rhs);
1333     changeInGuessed_ = rhs.changeInGuessed_;
1334   }
1335   return *this;
1336 }
1337 CbcBranchingObject *
clone() const1338 CbcIntegerPseudoCostBranchingObject::clone() const
1339 {
1340   return (new CbcIntegerPseudoCostBranchingObject(*this));
1341 }
1342 
1343 // Destructor
~CbcIntegerPseudoCostBranchingObject()1344 CbcIntegerPseudoCostBranchingObject::~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 */
1358 double
branch()1359 CbcIntegerPseudoCostBranchingObject::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 */
1373 CbcRangeCompare
compareBranchingObject(const CbcBranchingObject * brObj,const bool replaceIfOverlap)1374 CbcIntegerPseudoCostBranchingObject::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 */
CbcSwitchingBinary()1387 CbcSwitchingBinary::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 */
CbcSwitchingBinary(CbcSimpleIntegerDynamicPseudoCost * oldObject,int nOdd,const int * other,const int * otherRow)1401 CbcSwitchingBinary::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
CbcSwitchingBinary(const CbcSwitchingBinary & rhs)1508 CbcSwitchingBinary::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
1521 CbcObject *
clone() const1522 CbcSwitchingBinary::clone() const
1523 {
1524   return new CbcSwitchingBinary(*this);
1525 }
1526 
1527 // Assignment operator
1528 CbcSwitchingBinary &
operator =(const CbcSwitchingBinary & rhs)1529 CbcSwitchingBinary::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
~CbcSwitchingBinary()1547 CbcSwitchingBinary::~CbcSwitchingBinary()
1548 {
1549   delete[] zeroLowerBound_;
1550   delete[] otherVariable_;
1551 }
1552 // Add in zero switches
addZeroSwitches(int nAdd,const int * columns)1553 void 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)
same(const CbcSwitchingBinary * otherObject) const1580 bool CbcSwitchingBinary::same(const CbcSwitchingBinary *otherObject) const
1581 {
1582   bool okay = CbcSimpleIntegerDynamicPseudoCost::same(otherObject);
1583   return okay;
1584 }
1585 double
infeasibility(const OsiBranchingInformation * info,int & preferredWay) const1586 CbcSwitchingBinary::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
setAssociatedBounds(OsiSolverInterface * solver,int cleanBasis) const1629 int 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
checkAssociatedBounds(const OsiSolverInterface * solver,const double * solution,int printLevel,int state[3],int & nBadFixed) const1733 int 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 */
1815