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