1 | // Copyright (C) 2002, International Business Machines |
---|
2 | // Corporation and others. All Rights Reserved. |
---|
3 | #if defined(_MSC_VER) |
---|
4 | // Turn off compiler warning about long names |
---|
5 | # pragma warning(disable:4786) |
---|
6 | #endif |
---|
7 | #include <string> |
---|
8 | //#define CBC_DEBUG 1 |
---|
9 | //#define CHECK_CUT_COUNTS |
---|
10 | #include <cassert> |
---|
11 | #include <cfloat> |
---|
12 | #define CUTS |
---|
13 | #include "OsiSolverInterface.hpp" |
---|
14 | #include "CoinWarmStartBasis.hpp" |
---|
15 | #include "CoinTime.hpp" |
---|
16 | #include "CbcModel.hpp" |
---|
17 | #include "CbcNode.hpp" |
---|
18 | #include "CbcBranchActual.hpp" |
---|
19 | #include "OsiRowCut.hpp" |
---|
20 | #include "OsiRowCutDebugger.hpp" |
---|
21 | #include "OsiCuts.hpp" |
---|
22 | #include "CbcCountRowCut.hpp" |
---|
23 | #include "CbcMessage.hpp" |
---|
24 | #include "OsiClpSolverInterface.hpp" |
---|
25 | #include "ClpSimplexOther.hpp" |
---|
26 | using namespace std; |
---|
27 | #include "CglCutGenerator.hpp" |
---|
28 | // Default Constructor |
---|
29 | CbcNodeInfo::CbcNodeInfo () |
---|
30 | : |
---|
31 | numberPointingToThis_(0), |
---|
32 | parent_(NULL), |
---|
33 | owner_(NULL), |
---|
34 | numberCuts_(0), |
---|
35 | cuts_(NULL), |
---|
36 | numberRows_(0), |
---|
37 | numberBranchesLeft_(0) |
---|
38 | { |
---|
39 | #ifdef CHECK_NODE |
---|
40 | printf("CbcNodeInfo %x Constructor\n",this); |
---|
41 | #endif |
---|
42 | } |
---|
43 | // Constructor given parent |
---|
44 | CbcNodeInfo::CbcNodeInfo (CbcNodeInfo * parent) |
---|
45 | : |
---|
46 | numberPointingToThis_(2), |
---|
47 | parent_(parent), |
---|
48 | owner_(NULL), |
---|
49 | numberCuts_(0), |
---|
50 | cuts_(NULL), |
---|
51 | numberRows_(0), |
---|
52 | numberBranchesLeft_(2) |
---|
53 | { |
---|
54 | #ifdef CHECK_NODE |
---|
55 | printf("CbcNodeInfo %x Constructor from parent %x\n",this,parent_); |
---|
56 | #endif |
---|
57 | if (parent_) { |
---|
58 | numberRows_ = parent_->numberRows_+parent_->numberCuts_; |
---|
59 | } |
---|
60 | } |
---|
61 | // Copy Constructor |
---|
62 | CbcNodeInfo::CbcNodeInfo (const CbcNodeInfo & rhs) |
---|
63 | : |
---|
64 | numberPointingToThis_(rhs.numberPointingToThis_), |
---|
65 | parent_(rhs.parent_), |
---|
66 | owner_(rhs.owner_), |
---|
67 | numberCuts_(rhs.numberCuts_), |
---|
68 | cuts_(NULL), |
---|
69 | numberRows_(rhs.numberRows_), |
---|
70 | numberBranchesLeft_(rhs.numberBranchesLeft_) |
---|
71 | { |
---|
72 | #ifdef CHECK_NODE |
---|
73 | printf("CbcNodeInfo %x Copy constructor\n",this); |
---|
74 | #endif |
---|
75 | if (numberCuts_) { |
---|
76 | cuts_ = new CbcCountRowCut * [numberCuts_]; |
---|
77 | int n=0; |
---|
78 | for (int i=0;i<numberCuts_;i++) { |
---|
79 | CbcCountRowCut * thisCut = rhs.cuts_[i]; |
---|
80 | if (thisCut) { |
---|
81 | // I think this is correct - new one should take priority |
---|
82 | thisCut->setInfo(this,n); |
---|
83 | thisCut->increment(numberBranchesLeft_); |
---|
84 | cuts_[n++] = thisCut; |
---|
85 | } |
---|
86 | } |
---|
87 | numberCuts_=n; |
---|
88 | } |
---|
89 | } |
---|
90 | // Constructor given parent and owner |
---|
91 | CbcNodeInfo::CbcNodeInfo (CbcNodeInfo * parent, CbcNode * owner) |
---|
92 | : |
---|
93 | numberPointingToThis_(2), |
---|
94 | parent_(parent), |
---|
95 | owner_(owner), |
---|
96 | numberCuts_(0), |
---|
97 | cuts_(NULL), |
---|
98 | numberRows_(0), |
---|
99 | numberBranchesLeft_(2) |
---|
100 | { |
---|
101 | #ifdef CHECK_NODE |
---|
102 | printf("CbcNodeInfo %x Constructor from parent %x\n",this,parent_); |
---|
103 | #endif |
---|
104 | if (parent_) { |
---|
105 | numberRows_ = parent_->numberRows_+parent_->numberCuts_; |
---|
106 | } |
---|
107 | } |
---|
108 | |
---|
109 | /** |
---|
110 | Take care to detach from the owning CbcNode and decrement the reference |
---|
111 | count in the parent. If this is the last nodeInfo object pointing to the |
---|
112 | parent, make a recursive call to delete the parent. |
---|
113 | */ |
---|
114 | CbcNodeInfo::~CbcNodeInfo() |
---|
115 | { |
---|
116 | #ifdef CHECK_NODE |
---|
117 | printf("CbcNodeInfo %x Destructor parent %x\n",this,parent_); |
---|
118 | #endif |
---|
119 | |
---|
120 | assert(!numberPointingToThis_); |
---|
121 | // But they may be some left (max nodes?) |
---|
122 | for (int i=0;i<numberCuts_;i++) |
---|
123 | delete cuts_[i]; |
---|
124 | |
---|
125 | delete [] cuts_; |
---|
126 | if (owner_) |
---|
127 | owner_->nullNodeInfo(); |
---|
128 | if (parent_) { |
---|
129 | int numberLinks = parent_->decrement(); |
---|
130 | if (!numberLinks) delete parent_; |
---|
131 | } |
---|
132 | } |
---|
133 | |
---|
134 | |
---|
135 | //#define ALLCUTS |
---|
136 | void |
---|
137 | CbcNodeInfo::decrementCuts(int change) |
---|
138 | { |
---|
139 | int i; |
---|
140 | // get rid of all remaining if negative |
---|
141 | int changeThis; |
---|
142 | if (change<0) |
---|
143 | changeThis = numberBranchesLeft_; |
---|
144 | else |
---|
145 | changeThis = change; |
---|
146 | // decrement cut counts |
---|
147 | for (i=0;i<numberCuts_;i++) { |
---|
148 | if (cuts_[i]) { |
---|
149 | int number = cuts_[i]->decrement(changeThis); |
---|
150 | if (!number) { |
---|
151 | //printf("info %x del cut %d %x\n",this,i,cuts_[i]); |
---|
152 | delete cuts_[i]; |
---|
153 | cuts_[i]=NULL; |
---|
154 | } |
---|
155 | } |
---|
156 | } |
---|
157 | } |
---|
158 | void |
---|
159 | CbcNodeInfo::decrementParentCuts(int change) |
---|
160 | { |
---|
161 | if (parent_) { |
---|
162 | // get rid of all remaining if negative |
---|
163 | int changeThis; |
---|
164 | if (change<0) |
---|
165 | changeThis = numberBranchesLeft_; |
---|
166 | else |
---|
167 | changeThis = change; |
---|
168 | int i; |
---|
169 | // Get over-estimate of space needed for basis |
---|
170 | CoinWarmStartBasis dummy; |
---|
171 | dummy.setSize(0,numberRows_+numberCuts_); |
---|
172 | buildRowBasis(dummy); |
---|
173 | /* everything is zero (i.e. free) so we can use to see |
---|
174 | if latest basis */ |
---|
175 | CbcNodeInfo * thisInfo = parent_; |
---|
176 | while (thisInfo) |
---|
177 | thisInfo = thisInfo->buildRowBasis(dummy); |
---|
178 | // decrement cut counts |
---|
179 | thisInfo = parent_; |
---|
180 | int numberRows=numberRows_; |
---|
181 | while (thisInfo) { |
---|
182 | for (i=thisInfo->numberCuts_-1;i>=0;i--) { |
---|
183 | CoinWarmStartBasis::Status status = dummy.getArtifStatus(--numberRows); |
---|
184 | #ifdef ALLCUTS |
---|
185 | status = CoinWarmStartBasis::isFree; |
---|
186 | #endif |
---|
187 | if (thisInfo->cuts_[i]) { |
---|
188 | int number=1; |
---|
189 | if (status!=CoinWarmStartBasis::basic) { |
---|
190 | // tight - drop 1 or 2 |
---|
191 | if (change<0) |
---|
192 | number = thisInfo->cuts_[i]->decrement(changeThis); |
---|
193 | else |
---|
194 | number = thisInfo->cuts_[i]->decrement(change); |
---|
195 | } |
---|
196 | if (!number) { |
---|
197 | delete thisInfo->cuts_[i]; |
---|
198 | thisInfo->cuts_[i]=NULL; |
---|
199 | } |
---|
200 | } |
---|
201 | } |
---|
202 | thisInfo = thisInfo->parent_; |
---|
203 | } |
---|
204 | } |
---|
205 | } |
---|
206 | |
---|
207 | void |
---|
208 | CbcNodeInfo::incrementParentCuts(int change) |
---|
209 | { |
---|
210 | if (parent_) { |
---|
211 | int i; |
---|
212 | // Get over-estimate of space needed for basis |
---|
213 | CoinWarmStartBasis dummy; |
---|
214 | dummy.setSize(0,numberRows_+numberCuts_); |
---|
215 | /* everything is zero (i.e. free) so we can use to see |
---|
216 | if latest basis */ |
---|
217 | buildRowBasis(dummy); |
---|
218 | CbcNodeInfo * thisInfo = parent_; |
---|
219 | while (thisInfo) |
---|
220 | thisInfo = thisInfo->buildRowBasis(dummy); |
---|
221 | // increment cut counts |
---|
222 | thisInfo = parent_; |
---|
223 | int numberRows=numberRows_; |
---|
224 | while (thisInfo) { |
---|
225 | for (i=thisInfo->numberCuts_-1;i>=0;i--) { |
---|
226 | CoinWarmStartBasis::Status status = dummy.getArtifStatus(--numberRows); |
---|
227 | #ifdef ALLCUTS |
---|
228 | status = CoinWarmStartBasis::isFree; |
---|
229 | #endif |
---|
230 | if (thisInfo->cuts_[i]&&status!=CoinWarmStartBasis::basic) { |
---|
231 | thisInfo->cuts_[i]->increment(change); |
---|
232 | } |
---|
233 | } |
---|
234 | thisInfo = thisInfo->parent_; |
---|
235 | } |
---|
236 | } |
---|
237 | } |
---|
238 | |
---|
239 | /* |
---|
240 | Append cuts to the cuts_ array in a nodeInfo. The initial reference count |
---|
241 | is set to numberToBranchOn, which will normally be the number of arms |
---|
242 | defined for the CbcBranchingObject attached to the CbcNode that owns this |
---|
243 | CbcNodeInfo. |
---|
244 | */ |
---|
245 | void |
---|
246 | CbcNodeInfo::addCuts (OsiCuts & cuts, int numberToBranchOn, |
---|
247 | int * whichGenerator) |
---|
248 | { |
---|
249 | int numberCuts = cuts.sizeRowCuts(); |
---|
250 | if (numberCuts) { |
---|
251 | int i; |
---|
252 | if (!numberCuts_) { |
---|
253 | cuts_ = new CbcCountRowCut * [numberCuts]; |
---|
254 | } else { |
---|
255 | CbcCountRowCut ** temp = new CbcCountRowCut * [numberCuts+numberCuts_]; |
---|
256 | memcpy(temp,cuts_,numberCuts_*sizeof(CbcCountRowCut *)); |
---|
257 | delete [] cuts_; |
---|
258 | cuts_ = temp; |
---|
259 | } |
---|
260 | for (i=0;i<numberCuts;i++) { |
---|
261 | CbcCountRowCut * thisCut = new CbcCountRowCut(*cuts.rowCutPtr(i), |
---|
262 | this,numberCuts_); |
---|
263 | thisCut->increment(numberToBranchOn); |
---|
264 | cuts_[numberCuts_++] = thisCut; |
---|
265 | #ifdef CBC_DEBUG |
---|
266 | int n=thisCut->row().getNumElements(); |
---|
267 | #if CBC_DEBUG>1 |
---|
268 | printf("Cut %d has %d entries, rhs %g %g =>",i,n,thisCut->lb(), |
---|
269 | thisCut->ub()); |
---|
270 | #endif |
---|
271 | int j; |
---|
272 | #if CBC_DEBUG>1 |
---|
273 | const int * index = thisCut->row().getIndices(); |
---|
274 | #endif |
---|
275 | const double * element = thisCut->row().getElements(); |
---|
276 | for (j=0;j<n;j++) { |
---|
277 | #if CBC_DEBUG>1 |
---|
278 | printf(" (%d,%g)",index[j],element[j]); |
---|
279 | #endif |
---|
280 | assert(fabs(element[j])>1.00e-12); |
---|
281 | } |
---|
282 | printf("\n"); |
---|
283 | #endif |
---|
284 | } |
---|
285 | } |
---|
286 | } |
---|
287 | |
---|
288 | void |
---|
289 | CbcNodeInfo::addCuts(int numberCuts, CbcCountRowCut ** cut, |
---|
290 | int numberToBranchOn) |
---|
291 | { |
---|
292 | if (numberCuts) { |
---|
293 | int i; |
---|
294 | if (!numberCuts_) { |
---|
295 | cuts_ = new CbcCountRowCut * [numberCuts]; |
---|
296 | } else { |
---|
297 | CbcCountRowCut ** temp = new CbcCountRowCut * [numberCuts+numberCuts_]; |
---|
298 | memcpy(temp,cuts_,numberCuts_*sizeof(CbcCountRowCut *)); |
---|
299 | delete [] cuts_; |
---|
300 | cuts_ = temp; |
---|
301 | } |
---|
302 | for (i=0;i<numberCuts;i++) { |
---|
303 | CbcCountRowCut * thisCut = cut[i]; |
---|
304 | thisCut->setInfo(this,numberCuts_); |
---|
305 | //printf("info %x cut %d %x\n",this,i,thisCut); |
---|
306 | thisCut->increment(numberToBranchOn); |
---|
307 | cuts_[numberCuts_++] = thisCut; |
---|
308 | #ifdef CBC_DEBUG |
---|
309 | int n=thisCut->row().getNumElements(); |
---|
310 | #if CBC_DEBUG>1 |
---|
311 | printf("Cut %d has %d entries, rhs %g %g =>",i,n,thisCut->lb(), |
---|
312 | thisCut->ub()); |
---|
313 | #endif |
---|
314 | int j; |
---|
315 | #if CBC_DEBUG>1 |
---|
316 | const int * index = thisCut->row().getIndices(); |
---|
317 | #endif |
---|
318 | const double * element = thisCut->row().getElements(); |
---|
319 | for (j=0;j<n;j++) { |
---|
320 | #if CBC_DEBUG>1 |
---|
321 | printf(" (%d,%g)",index[j],element[j]); |
---|
322 | #endif |
---|
323 | assert(fabs(element[j])>1.00e-12); |
---|
324 | } |
---|
325 | printf("\n"); |
---|
326 | #endif |
---|
327 | } |
---|
328 | } |
---|
329 | } |
---|
330 | |
---|
331 | // delete cuts |
---|
332 | void |
---|
333 | CbcNodeInfo::deleteCuts(int numberToDelete, CbcCountRowCut ** cuts) |
---|
334 | { |
---|
335 | int i; |
---|
336 | int j; |
---|
337 | int last=-1; |
---|
338 | for (i=0;i<numberToDelete;i++) { |
---|
339 | CbcCountRowCut * next = cuts[i]; |
---|
340 | for (j=last+1;j<numberCuts_;j++) { |
---|
341 | if (next==cuts_[j]) |
---|
342 | break; |
---|
343 | } |
---|
344 | if (j==numberCuts_) { |
---|
345 | // start from beginning |
---|
346 | for (j=0;j<last;j++) { |
---|
347 | if (next==cuts_[j]) |
---|
348 | break; |
---|
349 | } |
---|
350 | assert(j<last); |
---|
351 | } |
---|
352 | last=j; |
---|
353 | int number = cuts_[j]->decrement(); |
---|
354 | if (!number) |
---|
355 | delete cuts_[j]; |
---|
356 | cuts_[j]=NULL; |
---|
357 | } |
---|
358 | j=0; |
---|
359 | for (i=0;i<numberCuts_;i++) { |
---|
360 | if (cuts_[i]) |
---|
361 | cuts_[j++]=cuts_[i]; |
---|
362 | } |
---|
363 | numberCuts_ = j; |
---|
364 | } |
---|
365 | |
---|
366 | // delete cuts |
---|
367 | void |
---|
368 | CbcNodeInfo::deleteCuts(int numberToDelete, int * which) |
---|
369 | { |
---|
370 | int i; |
---|
371 | for (i=0;i<numberToDelete;i++) { |
---|
372 | int iCut=which[i]; |
---|
373 | int number = cuts_[iCut]->decrement(); |
---|
374 | if (!number) |
---|
375 | delete cuts_[iCut]; |
---|
376 | cuts_[iCut]=NULL; |
---|
377 | } |
---|
378 | int j=0; |
---|
379 | for (i=0;i<numberCuts_;i++) { |
---|
380 | if (cuts_[i]) |
---|
381 | cuts_[j++]=cuts_[i]; |
---|
382 | } |
---|
383 | numberCuts_ = j; |
---|
384 | } |
---|
385 | |
---|
386 | // Really delete a cut |
---|
387 | void |
---|
388 | CbcNodeInfo::deleteCut(int whichOne) |
---|
389 | { |
---|
390 | assert(whichOne<numberCuts_); |
---|
391 | cuts_[whichOne]=NULL; |
---|
392 | } |
---|
393 | |
---|
394 | CbcFullNodeInfo::CbcFullNodeInfo() : |
---|
395 | CbcNodeInfo(), |
---|
396 | basis_(), |
---|
397 | numberIntegers_(0), |
---|
398 | lower_(NULL), |
---|
399 | upper_(NULL) |
---|
400 | { |
---|
401 | } |
---|
402 | CbcFullNodeInfo::CbcFullNodeInfo(CbcModel * model, |
---|
403 | int numberRowsAtContinuous) : |
---|
404 | CbcNodeInfo() |
---|
405 | { |
---|
406 | OsiSolverInterface * solver = model->solver(); |
---|
407 | numberRows_ = numberRowsAtContinuous; |
---|
408 | numberIntegers_ = model->numberIntegers(); |
---|
409 | int numberColumns = model->getNumCols(); |
---|
410 | lower_ = new double [numberColumns]; |
---|
411 | upper_ = new double [numberColumns]; |
---|
412 | const double * lower = solver->getColLower(); |
---|
413 | const double * upper = solver->getColUpper(); |
---|
414 | int i; |
---|
415 | |
---|
416 | for (i=0;i<numberColumns;i++) { |
---|
417 | lower_[i]=lower[i]; |
---|
418 | upper_[i]=upper[i]; |
---|
419 | } |
---|
420 | |
---|
421 | basis_ = dynamic_cast<CoinWarmStartBasis*>(solver->getWarmStart()); |
---|
422 | } |
---|
423 | |
---|
424 | CbcFullNodeInfo::CbcFullNodeInfo(const CbcFullNodeInfo & rhs) : |
---|
425 | CbcNodeInfo(rhs) |
---|
426 | { |
---|
427 | basis_= dynamic_cast<CoinWarmStartBasis *>(rhs.basis_->clone()) ; |
---|
428 | numberIntegers_=rhs.numberIntegers_; |
---|
429 | lower_=NULL; |
---|
430 | upper_=NULL; |
---|
431 | if (rhs.lower_!=NULL) { |
---|
432 | int numberColumns = basis_->getNumStructural(); |
---|
433 | lower_ = new double [numberColumns]; |
---|
434 | upper_ = new double [numberColumns]; |
---|
435 | assert (upper_!=NULL); |
---|
436 | memcpy(lower_,rhs.lower_,numberColumns*sizeof(double)); |
---|
437 | memcpy(upper_,rhs.upper_,numberColumns*sizeof(double)); |
---|
438 | } |
---|
439 | } |
---|
440 | |
---|
441 | CbcNodeInfo * |
---|
442 | CbcFullNodeInfo::clone() const |
---|
443 | { |
---|
444 | return (new CbcFullNodeInfo(*this)); |
---|
445 | } |
---|
446 | |
---|
447 | CbcFullNodeInfo::~CbcFullNodeInfo () |
---|
448 | { |
---|
449 | delete basis_ ; |
---|
450 | delete [] lower_; |
---|
451 | delete [] upper_; |
---|
452 | } |
---|
453 | |
---|
454 | /* |
---|
455 | The basis supplied as a parameter is deleted and replaced with a new basis |
---|
456 | appropriate for the node, and lower and upper bounds on variables are |
---|
457 | reset according to the stored bounds arrays. Any cuts associated with this |
---|
458 | node are added to the list in addCuts, but not actually added to the |
---|
459 | constraint system in the model. |
---|
460 | |
---|
461 | Why pass in a basis at all? The short answer is ``We need the parameter to |
---|
462 | pass out a basis, so might as well use it to pass in the size.'' |
---|
463 | |
---|
464 | A longer answer is that in practice we take a memory allocation hit up in |
---|
465 | addCuts1 (the only place applyToModel is called) when we setSize() the |
---|
466 | basis that's passed in. It's immediately tossed here in favour of a clone |
---|
467 | of the basis attached to this nodeInfo. This can probably be fixed, given |
---|
468 | a bit of thought. |
---|
469 | */ |
---|
470 | |
---|
471 | void CbcFullNodeInfo::applyToModel (CbcModel *model, |
---|
472 | CoinWarmStartBasis *&basis, |
---|
473 | CbcCountRowCut **addCuts, |
---|
474 | int ¤tNumberCuts) const |
---|
475 | |
---|
476 | { OsiSolverInterface *solver = model->solver() ; |
---|
477 | |
---|
478 | // branch - do bounds |
---|
479 | int i; |
---|
480 | int numberColumns = model->getNumCols(); |
---|
481 | for (i=0;i<numberColumns;i++) { |
---|
482 | solver->setColBounds(i,lower_[i],upper_[i]); |
---|
483 | } |
---|
484 | // move basis - but make sure size stays |
---|
485 | int numberRows = basis->getNumArtificial(); |
---|
486 | delete basis ; |
---|
487 | basis = dynamic_cast<CoinWarmStartBasis *>(basis_->clone()) ; |
---|
488 | basis->resize(numberRows,numberColumns); |
---|
489 | for (i=0;i<numberCuts_;i++) |
---|
490 | addCuts[currentNumberCuts+i]= cuts_[i]; |
---|
491 | currentNumberCuts += numberCuts_; |
---|
492 | assert(!parent_); |
---|
493 | return ; |
---|
494 | } |
---|
495 | |
---|
496 | /* Builds up row basis backwards (until original model). |
---|
497 | Returns NULL or previous one to apply . |
---|
498 | Depends on Free being 0 and impossible for cuts |
---|
499 | */ |
---|
500 | CbcNodeInfo * |
---|
501 | CbcFullNodeInfo::buildRowBasis(CoinWarmStartBasis & basis ) const |
---|
502 | { |
---|
503 | const unsigned int * saved = |
---|
504 | (const unsigned int *) basis_->getArtificialStatus(); |
---|
505 | unsigned int * now = |
---|
506 | (unsigned int *) basis.getArtificialStatus(); |
---|
507 | int number=basis_->getNumArtificial()>>4;; |
---|
508 | int i; |
---|
509 | for (i=0;i<number;i++) { |
---|
510 | if (!now[i]) |
---|
511 | now[i] = saved[i]; |
---|
512 | } |
---|
513 | return NULL; |
---|
514 | } |
---|
515 | |
---|
516 | // Default constructor |
---|
517 | CbcPartialNodeInfo::CbcPartialNodeInfo() |
---|
518 | |
---|
519 | : CbcNodeInfo(), |
---|
520 | basisDiff_(NULL), |
---|
521 | variables_(NULL), |
---|
522 | newBounds_(NULL), |
---|
523 | numberChangedBounds_(0) |
---|
524 | |
---|
525 | { /* this space intentionally left blank */ } |
---|
526 | |
---|
527 | // Constructor from current state |
---|
528 | CbcPartialNodeInfo::CbcPartialNodeInfo (CbcNodeInfo *parent, CbcNode *owner, |
---|
529 | int numberChangedBounds, |
---|
530 | const int *variables, |
---|
531 | const double *boundChanges, |
---|
532 | const CoinWarmStartDiff *basisDiff) |
---|
533 | : CbcNodeInfo(parent,owner) |
---|
534 | { |
---|
535 | basisDiff_ = basisDiff->clone() ; |
---|
536 | |
---|
537 | numberChangedBounds_ = numberChangedBounds; |
---|
538 | variables_ = new int [numberChangedBounds_]; |
---|
539 | newBounds_ = new double [numberChangedBounds_]; |
---|
540 | |
---|
541 | int i ; |
---|
542 | for (i=0;i<numberChangedBounds_;i++) { |
---|
543 | variables_[i]=variables[i]; |
---|
544 | newBounds_[i]=boundChanges[i]; |
---|
545 | } |
---|
546 | } |
---|
547 | |
---|
548 | CbcPartialNodeInfo::CbcPartialNodeInfo (const CbcPartialNodeInfo & rhs) |
---|
549 | |
---|
550 | : CbcNodeInfo(rhs.parent_) |
---|
551 | |
---|
552 | { basisDiff_ = rhs.basisDiff_->clone() ; |
---|
553 | |
---|
554 | numberChangedBounds_ = rhs.numberChangedBounds_; |
---|
555 | variables_ = new int [numberChangedBounds_]; |
---|
556 | newBounds_ = new double [numberChangedBounds_]; |
---|
557 | |
---|
558 | int i ; |
---|
559 | for (i=0;i<numberChangedBounds_;i++) { |
---|
560 | variables_[i]=rhs.variables_[i]; |
---|
561 | newBounds_[i]=rhs.newBounds_[i]; |
---|
562 | } |
---|
563 | } |
---|
564 | |
---|
565 | CbcNodeInfo * |
---|
566 | CbcPartialNodeInfo::clone() const |
---|
567 | { |
---|
568 | return (new CbcPartialNodeInfo(*this)); |
---|
569 | } |
---|
570 | |
---|
571 | |
---|
572 | CbcPartialNodeInfo::~CbcPartialNodeInfo () |
---|
573 | { |
---|
574 | delete basisDiff_ ; |
---|
575 | delete [] variables_; |
---|
576 | delete [] newBounds_; |
---|
577 | } |
---|
578 | |
---|
579 | |
---|
580 | /** |
---|
581 | The basis supplied as a parameter is incrementally modified, and lower and |
---|
582 | upper bounds on variables in the model are incrementally modified. Any |
---|
583 | cuts associated with this node are added to the list in addCuts. |
---|
584 | */ |
---|
585 | |
---|
586 | void CbcPartialNodeInfo::applyToModel (CbcModel *model, |
---|
587 | CoinWarmStartBasis *&basis, |
---|
588 | CbcCountRowCut **addCuts, |
---|
589 | int ¤tNumberCuts) const |
---|
590 | |
---|
591 | { OsiSolverInterface *solver = model->solver(); |
---|
592 | |
---|
593 | basis->applyDiff(basisDiff_) ; |
---|
594 | |
---|
595 | // branch - do bounds |
---|
596 | int i; |
---|
597 | for (i=0;i<numberChangedBounds_;i++) { |
---|
598 | int variable = variables_[i]; |
---|
599 | if ((variable&0x80000000)==0) { |
---|
600 | // lower bound changing |
---|
601 | solver->setColLower(variable,newBounds_[i]); |
---|
602 | } else { |
---|
603 | // upper bound changing |
---|
604 | solver->setColUpper(variable&0x7fffffff,newBounds_[i]); |
---|
605 | } |
---|
606 | } |
---|
607 | for (i=0;i<numberCuts_;i++) |
---|
608 | addCuts[currentNumberCuts+i]= cuts_[i]; |
---|
609 | currentNumberCuts += numberCuts_; |
---|
610 | return ; |
---|
611 | } |
---|
612 | |
---|
613 | /* Builds up row basis backwards (until original model). |
---|
614 | Returns NULL or previous one to apply . |
---|
615 | Depends on Free being 0 and impossible for cuts |
---|
616 | */ |
---|
617 | |
---|
618 | CbcNodeInfo * |
---|
619 | CbcPartialNodeInfo::buildRowBasis(CoinWarmStartBasis & basis ) const |
---|
620 | |
---|
621 | { basis.applyDiff(basisDiff_) ; |
---|
622 | |
---|
623 | return parent_ ; } |
---|
624 | |
---|
625 | |
---|
626 | CbcNode::CbcNode() : |
---|
627 | nodeInfo_(NULL), |
---|
628 | objectiveValue_(1.0e100), |
---|
629 | guessedObjectiveValue_(1.0e100), |
---|
630 | branch_(NULL), |
---|
631 | depth_(-1), |
---|
632 | numberUnsatisfied_(0), |
---|
633 | nodeNumber_(-1) |
---|
634 | { |
---|
635 | #ifdef CHECK_NODE |
---|
636 | printf("CbcNode %x Constructor\n",this); |
---|
637 | #endif |
---|
638 | } |
---|
639 | |
---|
640 | CbcNode::CbcNode(CbcModel * model, |
---|
641 | CbcNode * lastNode) : |
---|
642 | nodeInfo_(NULL), |
---|
643 | objectiveValue_(1.0e100), |
---|
644 | guessedObjectiveValue_(1.0e100), |
---|
645 | branch_(NULL), |
---|
646 | depth_(-1), |
---|
647 | numberUnsatisfied_(0), |
---|
648 | nodeNumber_(model->getNodeCount()) |
---|
649 | { |
---|
650 | #ifdef CHECK_NODE |
---|
651 | printf("CbcNode %x Constructor from model\n",this); |
---|
652 | #endif |
---|
653 | OsiSolverInterface * solver = model->solver(); |
---|
654 | objectiveValue_ = solver->getObjSense()*solver->getObjValue(); |
---|
655 | |
---|
656 | if (lastNode) |
---|
657 | lastNode->nodeInfo_->increment(); |
---|
658 | } |
---|
659 | |
---|
660 | |
---|
661 | void |
---|
662 | CbcNode::createInfo (CbcModel *model, |
---|
663 | CbcNode *lastNode, |
---|
664 | const CoinWarmStartBasis *lastws, |
---|
665 | const double *lastLower, const double *lastUpper, |
---|
666 | int numberOldActiveCuts,int numberNewCuts) |
---|
667 | { OsiSolverInterface * solver = model->solver(); |
---|
668 | /* |
---|
669 | The root --- no parent. Create full basis and bounds information. |
---|
670 | */ |
---|
671 | if (!lastNode) |
---|
672 | { nodeInfo_=new CbcFullNodeInfo(model,solver->getNumRows()); } |
---|
673 | /* |
---|
674 | Not the root. Create an edit from the parent's basis & bound information. |
---|
675 | This is not quite as straightforward as it seems. We need to reintroduce |
---|
676 | cuts we may have dropped out of the basis, in the correct position, because |
---|
677 | this whole process is strictly positional. Start by grabbing the current |
---|
678 | basis. |
---|
679 | */ |
---|
680 | else |
---|
681 | { const CoinWarmStartBasis* ws = |
---|
682 | dynamic_cast<const CoinWarmStartBasis*>(solver->getWarmStart()); |
---|
683 | assert(ws!=NULL); // make sure not volume |
---|
684 | //int numberArtificials = lastws->getNumArtificial(); |
---|
685 | int numberColumns = solver->getNumCols(); |
---|
686 | |
---|
687 | const double * lower = solver->getColLower(); |
---|
688 | const double * upper = solver->getColUpper(); |
---|
689 | |
---|
690 | int i; |
---|
691 | /* |
---|
692 | Create a clone and resize it to hold all the structural constraints, plus |
---|
693 | all the cuts: old cuts, both active and inactive (currentNumberCuts), and |
---|
694 | new cuts (numberNewCuts). |
---|
695 | |
---|
696 | TODO: You'd think that the set of constraints (logicals) in the expanded |
---|
697 | basis should match the set represented in lastws. At least, that's |
---|
698 | what I thought. But at the point I first looked hard at this bit of |
---|
699 | code, it turned out that lastws was the stripped basis produced at |
---|
700 | the end of addCuts(), rather than the raw basis handed back by |
---|
701 | addCuts1(). The expanded basis here is equivalent to the raw basis of |
---|
702 | addCuts1(). I said ``whoa, that's not good, I must have introduced a |
---|
703 | bug'' and went back to John's code to see where I'd gone wrong. |
---|
704 | And discovered the same `error' in his code. |
---|
705 | |
---|
706 | After a bit of thought, my conclusion is that correctness is not |
---|
707 | affected by whether lastws is the stripped or raw basis. The diffs |
---|
708 | have no semantics --- just a set of changes that need to be made |
---|
709 | to convert lastws into expanded. I think the only effect is that we |
---|
710 | store a lot more diffs (everything in expanded that's not covered by |
---|
711 | the stripped basis). But I need to give this more thought. There |
---|
712 | may well be some subtle error cases. |
---|
713 | |
---|
714 | In the mean time, I've twiddled addCuts() to set lastws to the raw |
---|
715 | basis. Makes me (Lou) less nervous to compare apples to apples. |
---|
716 | */ |
---|
717 | CoinWarmStartBasis *expanded = |
---|
718 | dynamic_cast<CoinWarmStartBasis *>(ws->clone()) ; |
---|
719 | int numberRowsAtContinuous = model->numberRowsAtContinuous(); |
---|
720 | int iFull = numberRowsAtContinuous+model->currentNumberCuts()+ |
---|
721 | numberNewCuts; |
---|
722 | //int numberArtificialsNow = iFull; |
---|
723 | //int maxBasisLength = ((iFull+15)>>4)+((numberColumns+15)>>4); |
---|
724 | //printf("l %d full %d\n",maxBasisLength,iFull); |
---|
725 | expanded->resize(iFull,numberColumns); |
---|
726 | #ifdef FULL_DEBUG |
---|
727 | printf("Before expansion: orig %d, old %d, new %d, current %d\n", |
---|
728 | numberRowsAtContinuous,numberOldActiveCuts,numberNewCuts, |
---|
729 | model->currentNumberCuts()) ; |
---|
730 | ws->print(); |
---|
731 | #endif |
---|
732 | /* |
---|
733 | Now fill in the expanded basis. Any indices beyond nPartial must |
---|
734 | be cuts created while processing this node --- they can be copied directly |
---|
735 | into the expanded basis. From nPartial down, pull the status of active cuts |
---|
736 | from ws, interleaving with a B entry for the deactivated (loose) cuts. |
---|
737 | */ |
---|
738 | int numberDropped = model->currentNumberCuts()-numberOldActiveCuts; |
---|
739 | int iCompact=iFull-numberDropped; |
---|
740 | CbcCountRowCut ** cut = model->addedCuts(); |
---|
741 | int nPartial = model->currentNumberCuts()+numberRowsAtContinuous; |
---|
742 | iFull--; |
---|
743 | for (;iFull>=nPartial;iFull--) { |
---|
744 | CoinWarmStartBasis::Status status = ws->getArtifStatus(--iCompact); |
---|
745 | //assert (status != CoinWarmStartBasis::basic); // may be permanent cut |
---|
746 | expanded->setArtifStatus(iFull,status); |
---|
747 | } |
---|
748 | for (;iFull>=numberRowsAtContinuous;iFull--) { |
---|
749 | if (cut[iFull-numberRowsAtContinuous]) { |
---|
750 | CoinWarmStartBasis::Status status = ws->getArtifStatus(--iCompact); |
---|
751 | // If no cut generator being used then we may have basic variables |
---|
752 | //if (model->getMaximumCutPasses()&& |
---|
753 | // status == CoinWarmStartBasis::basic) |
---|
754 | //printf("cut basic\n"); |
---|
755 | expanded->setArtifStatus(iFull,status); |
---|
756 | } else { |
---|
757 | expanded->setArtifStatus(iFull,CoinWarmStartBasis::basic); |
---|
758 | } |
---|
759 | } |
---|
760 | #ifdef FULL_DEBUG |
---|
761 | printf("Expanded basis\n"); |
---|
762 | expanded->print() ; |
---|
763 | printf("Diffing against\n") ; |
---|
764 | lastws->print() ; |
---|
765 | #endif |
---|
766 | /* |
---|
767 | Now that we have two bases in proper positional correspondence, creating |
---|
768 | the actual diff is dead easy. |
---|
769 | */ |
---|
770 | |
---|
771 | CoinWarmStartDiff *basisDiff = expanded->generateDiff(lastws) ; |
---|
772 | /* |
---|
773 | Diff the bound vectors. It's assumed the number of structural variables is |
---|
774 | not changing. Assuming that branching objects all involve integer variables, |
---|
775 | we should see at least one bound change as a consequence of processing this |
---|
776 | subproblem. Different types of branching objects could break this assertion. |
---|
777 | Yes - cuts will break |
---|
778 | */ |
---|
779 | double *boundChanges = new double [2*numberColumns] ; |
---|
780 | int *variables = new int [2*numberColumns] ; |
---|
781 | int numberChangedBounds=0; |
---|
782 | for (i=0;i<numberColumns;i++) { |
---|
783 | if (lower[i]!=lastLower[i]) { |
---|
784 | variables[numberChangedBounds]=i; |
---|
785 | boundChanges[numberChangedBounds++]=lower[i]; |
---|
786 | } |
---|
787 | if (upper[i]!=lastUpper[i]) { |
---|
788 | variables[numberChangedBounds]=i|0x80000000; |
---|
789 | boundChanges[numberChangedBounds++]=upper[i]; |
---|
790 | } |
---|
791 | #ifdef CBC_DEBUG |
---|
792 | if (lower[i]!=lastLower[i]) |
---|
793 | printf("lower on %d changed from %g to %g\n", |
---|
794 | i,lastLower[i],lower[i]); |
---|
795 | if (upper[i]!=lastUpper[i]) |
---|
796 | printf("upper on %d changed from %g to %g\n", |
---|
797 | i,lastUpper[i],upper[i]); |
---|
798 | #endif |
---|
799 | } |
---|
800 | #ifdef CBC_DEBUG |
---|
801 | printf("%d changed bounds\n",numberChangedBounds) ; |
---|
802 | #endif |
---|
803 | if (lastNode->branchingObject()->boundBranch()) |
---|
804 | assert (numberChangedBounds); |
---|
805 | /* |
---|
806 | Hand the lot over to the CbcPartialNodeInfo constructor, then clean up and |
---|
807 | return. |
---|
808 | */ |
---|
809 | nodeInfo_ = |
---|
810 | new CbcPartialNodeInfo(lastNode->nodeInfo_,this,numberChangedBounds, |
---|
811 | variables,boundChanges,basisDiff) ; |
---|
812 | delete basisDiff ; |
---|
813 | delete [] boundChanges; |
---|
814 | delete [] variables; |
---|
815 | delete expanded ; |
---|
816 | delete ws; |
---|
817 | } |
---|
818 | } |
---|
819 | |
---|
820 | /* |
---|
821 | The routine scans through the object list of the model looking for objects |
---|
822 | that indicate infeasibility. It tests each object using strong branching |
---|
823 | and selects the one with the least objective degradation. A corresponding |
---|
824 | branching object is left attached to lastNode. |
---|
825 | |
---|
826 | If strong branching is disabled, a candidate object is chosen essentially |
---|
827 | at random (whatever object ends up in pos'n 0 of the candidate array). |
---|
828 | |
---|
829 | If a branching candidate is found to be monotone, bounds are set to fix the |
---|
830 | variable and the routine immediately returns (the caller is expected to |
---|
831 | reoptimize). |
---|
832 | |
---|
833 | If a branching candidate is found to result in infeasibility in both |
---|
834 | directions, the routine immediately returns an indication of infeasibility. |
---|
835 | |
---|
836 | Returns: 0 both branch directions are feasible |
---|
837 | -1 branching variable is monotone |
---|
838 | -2 infeasible |
---|
839 | |
---|
840 | Original comments: |
---|
841 | Here could go cuts etc etc |
---|
842 | For now just fix on objective from strong branching. |
---|
843 | */ |
---|
844 | |
---|
845 | int CbcNode::chooseBranch (CbcModel *model, CbcNode *lastNode,int numberPassesLeft) |
---|
846 | |
---|
847 | { if (lastNode) |
---|
848 | depth_ = lastNode->depth_+1; |
---|
849 | else |
---|
850 | depth_ = 0; |
---|
851 | delete branch_; |
---|
852 | branch_=NULL; |
---|
853 | OsiSolverInterface * solver = model->solver(); |
---|
854 | double saveObjectiveValue = solver->getObjValue(); |
---|
855 | double objectiveValue = solver->getObjSense()*saveObjectiveValue; |
---|
856 | const double * lower = solver->getColLower(); |
---|
857 | const double * upper = solver->getColUpper(); |
---|
858 | int anyAction=0; |
---|
859 | double integerTolerance = |
---|
860 | model->getDblParam(CbcModel::CbcIntegerTolerance); |
---|
861 | int i; |
---|
862 | bool beforeSolution = model->getSolutionCount()==0; |
---|
863 | int numberStrong=model->numberStrong(); |
---|
864 | int saveNumberStrong=numberStrong; |
---|
865 | int numberObjects = model->numberObjects(); |
---|
866 | bool checkFeasibility = numberObjects>model->numberIntegers(); |
---|
867 | int maximumStrong = CoinMax(CoinMin(model->numberStrong(),numberObjects),1); |
---|
868 | int numberColumns = model->getNumCols(); |
---|
869 | double * saveUpper = new double[numberColumns]; |
---|
870 | double * saveLower = new double[numberColumns]; |
---|
871 | |
---|
872 | // Save solution in case heuristics need good solution later |
---|
873 | |
---|
874 | double * saveSolution = new double[numberColumns]; |
---|
875 | memcpy(saveSolution,solver->getColSolution(),numberColumns*sizeof(double)); |
---|
876 | |
---|
877 | /* |
---|
878 | Get a branching decision object. Use the default decision criteria unless |
---|
879 | the user has loaded a decision method into the model. |
---|
880 | */ |
---|
881 | CbcBranchDecision *decision = model->branchingMethod(); |
---|
882 | if (!decision) |
---|
883 | decision = new CbcBranchDefaultDecision(); |
---|
884 | |
---|
885 | typedef struct { |
---|
886 | CbcBranchingObject * possibleBranch; // what a branch would do |
---|
887 | double upMovement; // cost going up (and initial away from feasible) |
---|
888 | double downMovement; // cost going down |
---|
889 | int numIntInfeasUp ; // without odd ones |
---|
890 | int numObjInfeasUp ; // just odd ones |
---|
891 | bool finishedUp; // true if solver finished |
---|
892 | int numItersUp ; // number of iterations in solver |
---|
893 | int numIntInfeasDown ; // without odd ones |
---|
894 | int numObjInfeasDown ; // just odd ones |
---|
895 | bool finishedDown; // true if solver finished |
---|
896 | int numItersDown; // number of iterations in solver |
---|
897 | int objectNumber; // Which object it is |
---|
898 | int fix; // 0 if no fix, 1 if we can fix up, -1 if we can fix down |
---|
899 | } Strong; |
---|
900 | Strong * choice = new Strong[maximumStrong]; |
---|
901 | for (i=0;i<numberColumns;i++) { |
---|
902 | saveLower[i] = lower[i]; |
---|
903 | saveUpper[i] = upper[i]; |
---|
904 | } |
---|
905 | // May go round twice if strong branching fixes all local candidates |
---|
906 | bool finished=false; |
---|
907 | double estimatedDegradation=0.0; |
---|
908 | while(!finished) { |
---|
909 | finished=true; |
---|
910 | // Some objects may compute an estimate of best solution from here |
---|
911 | estimatedDegradation=0.0; |
---|
912 | int numberIntegerInfeasibilities=0; // without odd ones |
---|
913 | |
---|
914 | // We may go round this loop twice (only if we think we have solution) |
---|
915 | for (int iPass=0;iPass<2;iPass++) { |
---|
916 | |
---|
917 | // compute current state |
---|
918 | int numberObjectInfeasibilities; // just odd ones |
---|
919 | model->feasibleSolution( |
---|
920 | numberIntegerInfeasibilities, |
---|
921 | numberObjectInfeasibilities); |
---|
922 | // If forcePriority > 0 then we want best solution |
---|
923 | const double * bestSolution = NULL; |
---|
924 | int hotstartStrategy=model->getHotstartStrategy(); |
---|
925 | if (hotstartStrategy>0) { |
---|
926 | bestSolution = model->bestSolution(); |
---|
927 | } |
---|
928 | |
---|
929 | // Some objects may compute an estimate of best solution from here |
---|
930 | estimatedDegradation=0.0; |
---|
931 | numberUnsatisfied_ = 0; |
---|
932 | int bestPriority=INT_MAX; |
---|
933 | /* |
---|
934 | Scan for branching objects that indicate infeasibility. Choose the best |
---|
935 | maximumStrong candidates, using priority as the first criteria, then |
---|
936 | integer infeasibility. |
---|
937 | |
---|
938 | The algorithm is to fill the choice array with a set of good candidates (by |
---|
939 | infeasibility) with priority bestPriority. Finding a candidate with |
---|
940 | priority better (less) than bestPriority flushes the choice array. (This |
---|
941 | serves as initialization when the first candidate is found.) |
---|
942 | |
---|
943 | A new candidate is added to choices only if its infeasibility exceeds the |
---|
944 | current max infeasibility (mostAway). When a candidate is added, it |
---|
945 | replaces the candidate with the smallest infeasibility (tracked by |
---|
946 | iSmallest). |
---|
947 | */ |
---|
948 | int iSmallest = 0; |
---|
949 | double mostAway = 1.0e-100; |
---|
950 | for (i = 0 ; i < maximumStrong ; i++) choice[i].possibleBranch = NULL ; |
---|
951 | numberStrong=0; |
---|
952 | for (i=0;i<numberObjects;i++) { |
---|
953 | const CbcObject * object = model->object(i); |
---|
954 | int preferredWay; |
---|
955 | double infeasibility = object->infeasibility(preferredWay); |
---|
956 | int priorityLevel = object->priority(); |
---|
957 | if (bestSolution) { |
---|
958 | // we are doing hot start |
---|
959 | const CbcSimpleInteger * thisOne = dynamic_cast <const CbcSimpleInteger *> (object); |
---|
960 | if (thisOne) { |
---|
961 | int iColumn = thisOne->modelSequence(); |
---|
962 | if (saveUpper[iColumn]>saveLower[iColumn]) { |
---|
963 | double value = saveSolution[iColumn]; |
---|
964 | double targetValue = bestSolution[iColumn]; |
---|
965 | //double originalLower = thisOne->originalLower(); |
---|
966 | //double originalUpper = thisOne->originalUpper(); |
---|
967 | // switch off if not possible |
---|
968 | if (targetValue>=saveLower[iColumn]&&targetValue<=saveUpper[iColumn]) { |
---|
969 | /* priority outranks rest always if hotstartStrategy >1 |
---|
970 | otherwise can be downgraded if at correct level. |
---|
971 | Infeasibility may be increased by targetValue to choose 1.0 values first. |
---|
972 | */ |
---|
973 | if (fabs(value-targetValue)>integerTolerance) { |
---|
974 | if (value>targetValue) { |
---|
975 | infeasibility += value; |
---|
976 | preferredWay=-1; |
---|
977 | } else { |
---|
978 | infeasibility += targetValue; |
---|
979 | preferredWay=1; |
---|
980 | } |
---|
981 | } else if (hotstartStrategy>1) { |
---|
982 | if (targetValue==saveLower[iColumn]) { |
---|
983 | infeasibility += integerTolerance+1.0e-12; |
---|
984 | preferredWay=-1; |
---|
985 | } else if (targetValue==saveUpper[iColumn]) { |
---|
986 | infeasibility += integerTolerance+1.0e-12; |
---|
987 | preferredWay=1; |
---|
988 | } else { |
---|
989 | infeasibility += integerTolerance+1.0e-12; |
---|
990 | preferredWay=1; |
---|
991 | } |
---|
992 | } else { |
---|
993 | priorityLevel += 10000000; |
---|
994 | } |
---|
995 | } else { |
---|
996 | // switch off if not possible |
---|
997 | bestSolution=NULL; |
---|
998 | model->setHotstartStrategy(0); |
---|
999 | } |
---|
1000 | } |
---|
1001 | } |
---|
1002 | } |
---|
1003 | if (infeasibility) { |
---|
1004 | // Increase estimated degradation to solution |
---|
1005 | estimatedDegradation += CoinMin(object->upEstimate(),object->downEstimate()); |
---|
1006 | numberUnsatisfied_++; |
---|
1007 | // Better priority? Flush choices. |
---|
1008 | if (priorityLevel<bestPriority) { |
---|
1009 | int j; |
---|
1010 | iSmallest=0; |
---|
1011 | for (j=0;j<maximumStrong;j++) { |
---|
1012 | choice[j].upMovement=0.0; |
---|
1013 | delete choice[j].possibleBranch; |
---|
1014 | choice[j].possibleBranch=NULL; |
---|
1015 | } |
---|
1016 | bestPriority = priorityLevel; |
---|
1017 | mostAway=1.0e-100; |
---|
1018 | numberStrong=0; |
---|
1019 | } else if (priorityLevel>bestPriority) { |
---|
1020 | continue; |
---|
1021 | } |
---|
1022 | // Check for suitability based on infeasibility. |
---|
1023 | if (infeasibility>mostAway) { |
---|
1024 | //add to list |
---|
1025 | choice[iSmallest].upMovement=infeasibility; |
---|
1026 | delete choice[iSmallest].possibleBranch; |
---|
1027 | choice[iSmallest].possibleBranch=object->createBranch(preferredWay); |
---|
1028 | numberStrong = CoinMax(numberStrong,iSmallest+1); |
---|
1029 | // Save which object it was |
---|
1030 | choice[iSmallest].objectNumber=i; |
---|
1031 | int j; |
---|
1032 | iSmallest=-1; |
---|
1033 | mostAway = 1.0e50; |
---|
1034 | for (j=0;j<maximumStrong;j++) { |
---|
1035 | if (choice[j].upMovement<mostAway) { |
---|
1036 | mostAway=choice[j].upMovement; |
---|
1037 | iSmallest=j; |
---|
1038 | } |
---|
1039 | } |
---|
1040 | } |
---|
1041 | } |
---|
1042 | } |
---|
1043 | if (numberUnsatisfied_) { |
---|
1044 | // some infeasibilities - go to next steps |
---|
1045 | break; |
---|
1046 | } else if (!iPass) { |
---|
1047 | // looks like a solution - get paranoid |
---|
1048 | bool roundAgain=false; |
---|
1049 | // get basis |
---|
1050 | CoinWarmStartBasis * ws = dynamic_cast<CoinWarmStartBasis*>(solver->getWarmStart()); |
---|
1051 | if (!ws) |
---|
1052 | break; |
---|
1053 | for (i=0;i<numberColumns;i++) { |
---|
1054 | double value = saveSolution[i]; |
---|
1055 | if (value<lower[i]) { |
---|
1056 | saveSolution[i]=lower[i]; |
---|
1057 | roundAgain=true; |
---|
1058 | ws->setStructStatus(i,CoinWarmStartBasis::atLowerBound); |
---|
1059 | } else if (value>upper[i]) { |
---|
1060 | saveSolution[i]=upper[i]; |
---|
1061 | roundAgain=true; |
---|
1062 | ws->setStructStatus(i,CoinWarmStartBasis::atUpperBound); |
---|
1063 | } |
---|
1064 | } |
---|
1065 | if (roundAgain) { |
---|
1066 | // restore basis |
---|
1067 | solver->setWarmStart(ws); |
---|
1068 | delete ws; |
---|
1069 | solver->resolve(); |
---|
1070 | memcpy(saveSolution,solver->getColSolution(),numberColumns*sizeof(double)); |
---|
1071 | if (!solver->isProvenOptimal()) { |
---|
1072 | // infeasible |
---|
1073 | anyAction=-2; |
---|
1074 | break; |
---|
1075 | } |
---|
1076 | } else { |
---|
1077 | delete ws; |
---|
1078 | break; |
---|
1079 | } |
---|
1080 | } |
---|
1081 | } |
---|
1082 | /* Some solvers can do the strong branching calculations faster if |
---|
1083 | they do them all at once. At present only Clp does for ordinary |
---|
1084 | integers but I think this coding would be easy to modify |
---|
1085 | */ |
---|
1086 | bool allNormal=true; // to say if we can do fast strong branching |
---|
1087 | // Say which one will be best |
---|
1088 | int bestChoice=0; |
---|
1089 | double worstInfeasibility=0.0; |
---|
1090 | for (i=0;i<numberStrong;i++) { |
---|
1091 | choice[i].numIntInfeasUp = numberUnsatisfied_; |
---|
1092 | choice[i].numIntInfeasDown = numberUnsatisfied_; |
---|
1093 | choice[i].fix=0; // say not fixed |
---|
1094 | if (!dynamic_cast <const CbcSimpleInteger *> (model->object(choice[i].objectNumber))) |
---|
1095 | allNormal=false; // Something odd so lets skip clever fast branching |
---|
1096 | if ( !model->object(choice[i].objectNumber)->boundBranch()) |
---|
1097 | numberStrong=0; // switch off |
---|
1098 | // Do best choice in case switched off |
---|
1099 | if (choice[i].upMovement>worstInfeasibility) { |
---|
1100 | worstInfeasibility=choice[i].upMovement; |
---|
1101 | bestChoice=i; |
---|
1102 | } |
---|
1103 | } |
---|
1104 | // If we have hit max time don't do strong branching |
---|
1105 | bool hitMaxTime = ( CoinCpuTime()-model->getDblParam(CbcModel::CbcStartSeconds) > |
---|
1106 | model->getDblParam(CbcModel::CbcMaximumSeconds)); |
---|
1107 | // also give up if we are looping round too much |
---|
1108 | if (hitMaxTime||numberPassesLeft<=0) |
---|
1109 | numberStrong=0; |
---|
1110 | /* |
---|
1111 | Is strong branching enabled? If so, set up and do it. Otherwise, we'll |
---|
1112 | fall through to simple branching. |
---|
1113 | |
---|
1114 | Setup for strong branching involves saving the current basis (for restoration |
---|
1115 | afterwards) and setting up for hot starts. |
---|
1116 | */ |
---|
1117 | if (numberStrong&&saveNumberStrong) { |
---|
1118 | |
---|
1119 | bool solveAll=false; // set true to say look at all even if some fixed (experiment) |
---|
1120 | solveAll=true; |
---|
1121 | // worth trying if too many times |
---|
1122 | // Save basis |
---|
1123 | CoinWarmStart * ws = solver->getWarmStart(); |
---|
1124 | // save limit |
---|
1125 | int saveLimit; |
---|
1126 | solver->getIntParam(OsiMaxNumIterationHotStart,saveLimit); |
---|
1127 | if (beforeSolution&&saveLimit<100) |
---|
1128 | solver->setIntParam(OsiMaxNumIterationHotStart,100); // go to end |
---|
1129 | |
---|
1130 | /* If we are doing all strong branching in one go then we create new arrays |
---|
1131 | to store information. If clp NULL then doing old way. |
---|
1132 | Going down - |
---|
1133 | outputSolution[2*i] is final solution. |
---|
1134 | outputStuff[2*i] is status (0 - finished, 1 infeas, other unknown |
---|
1135 | outputStuff[2*i+numberStrong] is number iterations |
---|
1136 | On entry newUpper[i] is new upper bound, on exit obj change |
---|
1137 | Going up - |
---|
1138 | outputSolution[2*i+1] is final solution. |
---|
1139 | outputStuff[2*i+1] is status (0 - finished, 1 infeas, other unknown |
---|
1140 | outputStuff[2*i+1+numberStrong] is number iterations |
---|
1141 | On entry newLower[i] is new lower bound, on exit obj change |
---|
1142 | */ |
---|
1143 | OsiClpSolverInterface * osiclp = dynamic_cast< OsiClpSolverInterface*> (solver); |
---|
1144 | ClpSimplex * clp=NULL; |
---|
1145 | double * newLower = NULL; |
---|
1146 | double * newUpper = NULL; |
---|
1147 | double ** outputSolution=NULL; |
---|
1148 | int * outputStuff=NULL; |
---|
1149 | // Go back to normal way if user wants it |
---|
1150 | if (osiclp&&(osiclp->specialOptions()&16)!=0&&osiclp->specialOptions()>0) |
---|
1151 | allNormal=false; |
---|
1152 | if (osiclp&& allNormal) { |
---|
1153 | clp = osiclp->getModelPtr(); |
---|
1154 | // Clp - do a different way |
---|
1155 | newLower = new double[numberStrong]; |
---|
1156 | newUpper = new double[numberStrong]; |
---|
1157 | outputSolution = new double * [2*numberStrong]; |
---|
1158 | outputStuff = new int [4*numberStrong]; |
---|
1159 | int * which = new int[numberStrong]; |
---|
1160 | int startFinishOptions; |
---|
1161 | int specialOptions = osiclp->specialOptions(); |
---|
1162 | int clpOptions = clp->specialOptions(); |
---|
1163 | int returnCode=0; |
---|
1164 | #define CRUNCH |
---|
1165 | #ifdef CRUNCH |
---|
1166 | // Crunch down problem |
---|
1167 | int numberRows = clp->numberRows(); |
---|
1168 | // Use dual region |
---|
1169 | double * rhs = clp->dualRowSolution(); |
---|
1170 | int * whichRow = new int[3*numberRows]; |
---|
1171 | int * whichColumn = new int[2*numberColumns]; |
---|
1172 | int nBound; |
---|
1173 | ClpSimplex * small = ((ClpSimplexOther *) clp)->crunch(rhs,whichRow,whichColumn,nBound,true); |
---|
1174 | if (!small) { |
---|
1175 | anyAction=-2; |
---|
1176 | //printf("XXXX Inf by inspection\n"); |
---|
1177 | delete [] whichColumn; |
---|
1178 | whichColumn=NULL; |
---|
1179 | delete [] whichRow; |
---|
1180 | whichRow=NULL; |
---|
1181 | break; |
---|
1182 | } else { |
---|
1183 | clp = small; |
---|
1184 | } |
---|
1185 | #else |
---|
1186 | int saveLogLevel = clp->logLevel(); |
---|
1187 | int saveMaxIts = clp->maximumIterations(); |
---|
1188 | #endif |
---|
1189 | clp->setLogLevel(0); |
---|
1190 | if((specialOptions&1)==0) { |
---|
1191 | startFinishOptions=0; |
---|
1192 | clp->setSpecialOptions(clpOptions|(64|1024)); |
---|
1193 | } else { |
---|
1194 | startFinishOptions=1+2+4; |
---|
1195 | //startFinishOptions=1+4; // for moment re-factorize |
---|
1196 | if((specialOptions&4)==0) |
---|
1197 | clp->setSpecialOptions(clpOptions|(64|128|512|1024|4096)); |
---|
1198 | else |
---|
1199 | clp->setSpecialOptions(clpOptions|(64|128|512|1024|2048|4096)); |
---|
1200 | } |
---|
1201 | // User may want to clean up before strong branching |
---|
1202 | if ((clp->specialOptions()&32)!=0) { |
---|
1203 | clp->primal(1); |
---|
1204 | if (clp->numberIterations()) |
---|
1205 | model->messageHandler()->message(CBC_ITERATE_STRONG,model->messages()) |
---|
1206 | << clp->numberIterations() |
---|
1207 | <<CoinMessageEol; |
---|
1208 | } |
---|
1209 | clp->setMaximumIterations(saveLimit); |
---|
1210 | #ifdef CRUNCH |
---|
1211 | int * backColumn = whichColumn+numberColumns; |
---|
1212 | #endif |
---|
1213 | for (i=0;i<numberStrong;i++) { |
---|
1214 | int iObject = choice[i].objectNumber; |
---|
1215 | const CbcObject * object = model->object(iObject); |
---|
1216 | const CbcSimpleInteger * simple = dynamic_cast <const CbcSimpleInteger *> (object); |
---|
1217 | int iSequence = simple->modelSequence(); |
---|
1218 | newLower[i]= ceil(saveSolution[iSequence]); |
---|
1219 | newUpper[i]= floor(saveSolution[iSequence]); |
---|
1220 | #ifdef CRUNCH |
---|
1221 | iSequence = backColumn[iSequence]; |
---|
1222 | assert (iSequence>=0); |
---|
1223 | #endif |
---|
1224 | which[i]=iSequence; |
---|
1225 | outputSolution[2*i]= new double [numberColumns]; |
---|
1226 | outputSolution[2*i+1]= new double [numberColumns]; |
---|
1227 | } |
---|
1228 | //clp->writeMps("bad"); |
---|
1229 | returnCode=clp->strongBranching(numberStrong,which, |
---|
1230 | newLower, newUpper,outputSolution, |
---|
1231 | outputStuff,outputStuff+2*numberStrong,!solveAll,false, |
---|
1232 | startFinishOptions); |
---|
1233 | #ifndef CRUNCH |
---|
1234 | clp->setSpecialOptions(clpOptions); // restore |
---|
1235 | clp->setMaximumIterations(saveMaxIts); |
---|
1236 | clp->setLogLevel(saveLogLevel); |
---|
1237 | #endif |
---|
1238 | if (returnCode==-2) { |
---|
1239 | // bad factorization!!! |
---|
1240 | // Doing normal way |
---|
1241 | // Mark hot start |
---|
1242 | solver->markHotStart(); |
---|
1243 | clp = NULL; |
---|
1244 | } else { |
---|
1245 | #ifdef CRUNCH |
---|
1246 | // extract solution |
---|
1247 | //bool checkSol=true; |
---|
1248 | for (i=0;i<numberStrong;i++) { |
---|
1249 | int iObject = choice[i].objectNumber; |
---|
1250 | const CbcObject * object = model->object(iObject); |
---|
1251 | const CbcSimpleInteger * simple = dynamic_cast <const CbcSimpleInteger *> (object); |
---|
1252 | int iSequence = simple->modelSequence(); |
---|
1253 | which[i]=iSequence; |
---|
1254 | double * sol = outputSolution[2*i]; |
---|
1255 | double * sol2 = outputSolution[2*i+1]; |
---|
1256 | //bool x=true; |
---|
1257 | //bool x2=true; |
---|
1258 | for (int iColumn=numberColumns-1;iColumn>=0;iColumn--) { |
---|
1259 | int jColumn = backColumn[iColumn]; |
---|
1260 | if (jColumn>=0) { |
---|
1261 | sol[iColumn]=sol[jColumn]; |
---|
1262 | sol2[iColumn]=sol2[jColumn]; |
---|
1263 | } else { |
---|
1264 | sol[iColumn]=saveSolution[iColumn]; |
---|
1265 | sol2[iColumn]=saveSolution[iColumn]; |
---|
1266 | } |
---|
1267 | //printf("strong %d col %d sol %g sol2 %g\n",i,iColumn,sol[iColumn], |
---|
1268 | // sol2[iColumn]); |
---|
1269 | //if (fabs(sol[iColumn]-0.5)<0.499) |
---|
1270 | //x=false; |
---|
1271 | //if (fabs(sol2[iColumn]-0.5)<0.499) |
---|
1272 | //x2=false; |
---|
1273 | } |
---|
1274 | #if 0 |
---|
1275 | if( outputStuff[2*i]!=0&&outputStuff[2*i+1]!=0) |
---|
1276 | checkSol=false; |
---|
1277 | if (x&&checkSol) { |
---|
1278 | int iStatus = outputStuff[2*i]; |
---|
1279 | double newObjectiveValue = objectiveValue+newUpper[i]; |
---|
1280 | //printf("solution found - status %d obj %g\n",iStatus,newObjectiveValue); |
---|
1281 | if (!iStatus) { |
---|
1282 | const double * lower = solver->getRowLower(); |
---|
1283 | const double * upper = solver->getRowUpper(); |
---|
1284 | int numberRows = solver->getNumRows(); |
---|
1285 | double * rhs = new double [numberRows]; |
---|
1286 | CoinZeroN(rhs,numberRows); |
---|
1287 | { |
---|
1288 | int numberColumns = solver->getNumCols(); |
---|
1289 | const double * columnLower = solver->getColLower(); |
---|
1290 | const double * columnUpper = solver->getColUpper(); |
---|
1291 | int numberRows = solver->getNumRows(); |
---|
1292 | |
---|
1293 | const double * element = matrix->getElements(); |
---|
1294 | const int * row = matrix->getIndices(); |
---|
1295 | const CoinBigIndex * columnStart = matrix->getVectorStarts(); |
---|
1296 | const int * columnLength = matrix->getVectorLengths(); |
---|
1297 | CoinZeroN(rhs,numberRows); |
---|
1298 | int iColumn; |
---|
1299 | for (iColumn=0;iColumn<numberColumns;iColumn++) { |
---|
1300 | double lower = columnLower[iColumn]; |
---|
1301 | double upper = columnUpper[iColumn]; |
---|
1302 | double solValue = sol[iColumn]; |
---|
1303 | assert (solValue>=lower-1.0e-4&&solValue<upper+1.0e-4); |
---|
1304 | for (CoinBigIndex j = columnStart[iColumn]; |
---|
1305 | j<columnStart[iColumn]+columnLength[iColumn];j++) { |
---|
1306 | int iRow = row[j]; |
---|
1307 | double value = element[j]; |
---|
1308 | rhs[iRow] += solValue*value; |
---|
1309 | if (iRow==-19) |
---|
1310 | printf("col %d has sol %g and el %g , rhs now %g\n", |
---|
1311 | iColumn,solValue,element[j],rhs[19]); |
---|
1312 | } |
---|
1313 | } |
---|
1314 | } |
---|
1315 | for (int i=0;i<numberRows;i++) { |
---|
1316 | assert (rhs[i]>lower[i]-1.0e-3); |
---|
1317 | assert (rhs[i]<upper[i]+1.0e-3); |
---|
1318 | } |
---|
1319 | delete [] rhs; |
---|
1320 | } |
---|
1321 | } |
---|
1322 | if (x2&&checkSol) { |
---|
1323 | int iStatus = outputStuff[2*i+1]; |
---|
1324 | double newObjectiveValue = objectiveValue+newLower[i]; |
---|
1325 | //printf("solution found - status %d obj %g\n",iStatus,newObjectiveValue); |
---|
1326 | if (!iStatus) { |
---|
1327 | const double * lower = solver->getRowLower(); |
---|
1328 | const double * upper = solver->getRowUpper(); |
---|
1329 | int numberRows = solver->getNumRows(); |
---|
1330 | double * rhs = new double [numberRows]; |
---|
1331 | CoinZeroN(rhs,numberRows); |
---|
1332 | solver->getMatrixByCol()->times(sol2,rhs) ; |
---|
1333 | for (int i=0;i<numberRows;i++) { |
---|
1334 | assert (rhs[i]>lower[i]-1.0e-4); |
---|
1335 | assert (rhs[i]<upper[i]+1.0e-4); |
---|
1336 | } |
---|
1337 | delete [] rhs; |
---|
1338 | } |
---|
1339 | } |
---|
1340 | #endif |
---|
1341 | } |
---|
1342 | #endif |
---|
1343 | } |
---|
1344 | #ifdef CRUNCH |
---|
1345 | delete [] whichColumn; |
---|
1346 | delete [] whichRow; |
---|
1347 | delete small; |
---|
1348 | #endif |
---|
1349 | delete [] which; |
---|
1350 | } else { |
---|
1351 | // Doing normal way |
---|
1352 | // Mark hot start |
---|
1353 | solver->markHotStart(); |
---|
1354 | } |
---|
1355 | /* |
---|
1356 | Open a loop to do the strong branching LPs. For each candidate variable, |
---|
1357 | solve an LP with the variable forced down, then up. If a direction turns |
---|
1358 | out to be infeasible or monotonic (i.e., over the dual objective cutoff), |
---|
1359 | force the objective change to be big (1.0e100). If we determine the problem |
---|
1360 | is infeasible, or find a monotone variable, escape the loop. |
---|
1361 | |
---|
1362 | TODO: The `restore bounds' part might be better encapsulated as an |
---|
1363 | unbranch() method. Branching objects more exotic than simple integers |
---|
1364 | or cliques might not restrict themselves to variable bounds. |
---|
1365 | |
---|
1366 | TODO: Virtuous solvers invalidate the current solution (or give bogus |
---|
1367 | results :-) when the bounds are changed out from under them. So we |
---|
1368 | need to do all the work associated with finding a new solution before |
---|
1369 | restoring the bounds. |
---|
1370 | */ |
---|
1371 | for (i = 0 ; i < numberStrong ; i++) |
---|
1372 | { double objectiveChange ; |
---|
1373 | double newObjectiveValue=1.0e100; |
---|
1374 | // status is 0 finished, 1 infeasible and other |
---|
1375 | int iStatus; |
---|
1376 | /* |
---|
1377 | Try the down direction first. (Specify the initial branching alternative as |
---|
1378 | down with a call to way(-1). Each subsequent call to branch() performs the |
---|
1379 | specified branch and advances the branch object state to the next branch |
---|
1380 | alternative.) |
---|
1381 | */ |
---|
1382 | if (!clp) { |
---|
1383 | choice[i].possibleBranch->way(-1) ; |
---|
1384 | choice[i].possibleBranch->branch() ; |
---|
1385 | bool feasible=true; |
---|
1386 | if (checkFeasibility) { |
---|
1387 | // check branching did not make infeasible |
---|
1388 | int iColumn; |
---|
1389 | int numberColumns = solver->getNumCols(); |
---|
1390 | const double * columnLower = solver->getColLower(); |
---|
1391 | const double * columnUpper = solver->getColUpper(); |
---|
1392 | for (iColumn= 0;iColumn<numberColumns;iColumn++) { |
---|
1393 | if (columnLower[iColumn]>columnUpper[iColumn]+1.0e-5) |
---|
1394 | feasible=false; |
---|
1395 | } |
---|
1396 | } |
---|
1397 | if (feasible) { |
---|
1398 | solver->solveFromHotStart() ; |
---|
1399 | /* |
---|
1400 | We now have an estimate of objective degradation that we can use for strong |
---|
1401 | branching. If we're over the cutoff, the variable is monotone up. |
---|
1402 | If we actually made it to optimality, check for a solution, and if we have |
---|
1403 | a good one, call setBestSolution to process it. Note that this may reduce the |
---|
1404 | cutoff, so we check again to see if we can declare this variable monotone. |
---|
1405 | */ |
---|
1406 | if (solver->isProvenOptimal()) |
---|
1407 | iStatus=0; // optimal |
---|
1408 | else if (solver->isIterationLimitReached() |
---|
1409 | &&!solver->isDualObjectiveLimitReached()) |
---|
1410 | iStatus=2; // unknown |
---|
1411 | else |
---|
1412 | iStatus=1; // infeasible |
---|
1413 | newObjectiveValue = solver->getObjSense()*solver->getObjValue(); |
---|
1414 | choice[i].numItersDown = solver->getIterationCount(); |
---|
1415 | } else { |
---|
1416 | iStatus=1; // infeasible |
---|
1417 | newObjectiveValue = 1.0e100; |
---|
1418 | choice[i].numItersDown = 0; |
---|
1419 | } |
---|
1420 | objectiveChange = newObjectiveValue-objectiveValue ; |
---|
1421 | } else { |
---|
1422 | iStatus = outputStuff[2*i]; |
---|
1423 | choice[i].numItersDown = outputStuff[2*numberStrong+2*i]; |
---|
1424 | newObjectiveValue = objectiveValue+newUpper[i]; |
---|
1425 | solver->setColSolution(outputSolution[2*i]); |
---|
1426 | } |
---|
1427 | objectiveChange = newObjectiveValue - objectiveValue; |
---|
1428 | if (!iStatus) { |
---|
1429 | choice[i].finishedDown = true ; |
---|
1430 | if (newObjectiveValue>=model->getCutoff()) { |
---|
1431 | objectiveChange = 1.0e100; // say infeasible |
---|
1432 | } else { |
---|
1433 | // See if integer solution |
---|
1434 | if (model->feasibleSolution(choice[i].numIntInfeasDown, |
---|
1435 | choice[i].numObjInfeasDown)) |
---|
1436 | { model->setBestSolution(CBC_STRONGSOL, |
---|
1437 | newObjectiveValue, |
---|
1438 | solver->getColSolution()) ; |
---|
1439 | if (newObjectiveValue >= model->getCutoff()) // *new* cutoff |
---|
1440 | objectiveChange = 1.0e100 ; |
---|
1441 | } |
---|
1442 | } |
---|
1443 | } else if (iStatus==1) { |
---|
1444 | objectiveChange = 1.0e100 ; |
---|
1445 | } else { |
---|
1446 | // Can't say much as we did not finish |
---|
1447 | choice[i].finishedDown = false ; |
---|
1448 | } |
---|
1449 | choice[i].downMovement = objectiveChange ; |
---|
1450 | |
---|
1451 | // restore bounds |
---|
1452 | if (!clp) |
---|
1453 | { for (int j=0;j<numberColumns;j++) { |
---|
1454 | if (saveLower[j] != lower[j]) |
---|
1455 | solver->setColLower(j,saveLower[j]); |
---|
1456 | if (saveUpper[j] != upper[j]) |
---|
1457 | solver->setColUpper(j,saveUpper[j]); |
---|
1458 | } |
---|
1459 | } |
---|
1460 | //printf("Down on %d, status is %d, obj %g its %d cost %g finished %d inf %d infobj %d\n", |
---|
1461 | // choice[i].objectNumber,iStatus,newObjectiveValue,choice[i].numItersDown, |
---|
1462 | // choice[i].downMovement,choice[i].finishedDown,choice[i].numIntInfeasDown, |
---|
1463 | // choice[i].numObjInfeasDown); |
---|
1464 | |
---|
1465 | // repeat the whole exercise, forcing the variable up |
---|
1466 | if (!clp) { |
---|
1467 | choice[i].possibleBranch->branch(); |
---|
1468 | bool feasible=true; |
---|
1469 | if (checkFeasibility) { |
---|
1470 | // check branching did not make infeasible |
---|
1471 | int iColumn; |
---|
1472 | int numberColumns = solver->getNumCols(); |
---|
1473 | const double * columnLower = solver->getColLower(); |
---|
1474 | const double * columnUpper = solver->getColUpper(); |
---|
1475 | for (iColumn= 0;iColumn<numberColumns;iColumn++) { |
---|
1476 | if (columnLower[iColumn]>columnUpper[iColumn]+1.0e-5) |
---|
1477 | feasible=false; |
---|
1478 | } |
---|
1479 | } |
---|
1480 | if (feasible) { |
---|
1481 | solver->solveFromHotStart() ; |
---|
1482 | /* |
---|
1483 | We now have an estimate of objective degradation that we can use for strong |
---|
1484 | branching. If we're over the cutoff, the variable is monotone up. |
---|
1485 | If we actually made it to optimality, check for a solution, and if we have |
---|
1486 | a good one, call setBestSolution to process it. Note that this may reduce the |
---|
1487 | cutoff, so we check again to see if we can declare this variable monotone. |
---|
1488 | */ |
---|
1489 | if (solver->isProvenOptimal()) |
---|
1490 | iStatus=0; // optimal |
---|
1491 | else if (solver->isIterationLimitReached() |
---|
1492 | &&!solver->isDualObjectiveLimitReached()) |
---|
1493 | iStatus=2; // unknown |
---|
1494 | else |
---|
1495 | iStatus=1; // infeasible |
---|
1496 | newObjectiveValue = solver->getObjSense()*solver->getObjValue(); |
---|
1497 | choice[i].numItersUp = solver->getIterationCount(); |
---|
1498 | } else { |
---|
1499 | iStatus=1; // infeasible |
---|
1500 | newObjectiveValue = 1.0e100; |
---|
1501 | choice[i].numItersDown = 0; |
---|
1502 | } |
---|
1503 | objectiveChange = newObjectiveValue-objectiveValue ; |
---|
1504 | } else { |
---|
1505 | iStatus = outputStuff[2*i+1]; |
---|
1506 | choice[i].numItersUp = outputStuff[2*numberStrong+2*i+1]; |
---|
1507 | newObjectiveValue = objectiveValue+newLower[i]; |
---|
1508 | solver->setColSolution(outputSolution[2*i+1]); |
---|
1509 | } |
---|
1510 | objectiveChange = newObjectiveValue - objectiveValue; |
---|
1511 | if (!iStatus) { |
---|
1512 | choice[i].finishedUp = true ; |
---|
1513 | if (newObjectiveValue>=model->getCutoff()) { |
---|
1514 | objectiveChange = 1.0e100; // say infeasible |
---|
1515 | } else { |
---|
1516 | // See if integer solution |
---|
1517 | if (model->feasibleSolution(choice[i].numIntInfeasUp, |
---|
1518 | choice[i].numObjInfeasUp)) |
---|
1519 | { model->setBestSolution(CBC_STRONGSOL, |
---|
1520 | newObjectiveValue, |
---|
1521 | solver->getColSolution()) ; |
---|
1522 | if (newObjectiveValue >= model->getCutoff()) // *new* cutoff |
---|
1523 | objectiveChange = 1.0e100 ; |
---|
1524 | } |
---|
1525 | } |
---|
1526 | } else if (iStatus==1) { |
---|
1527 | objectiveChange = 1.0e100 ; |
---|
1528 | } else { |
---|
1529 | // Can't say much as we did not finish |
---|
1530 | choice[i].finishedUp = false ; |
---|
1531 | } |
---|
1532 | choice[i].upMovement = objectiveChange ; |
---|
1533 | |
---|
1534 | // restore bounds |
---|
1535 | if (!clp) |
---|
1536 | { for (int j=0;j<numberColumns;j++) { |
---|
1537 | if (saveLower[j] != lower[j]) |
---|
1538 | solver->setColLower(j,saveLower[j]); |
---|
1539 | if (saveUpper[j] != upper[j]) |
---|
1540 | solver->setColUpper(j,saveUpper[j]); |
---|
1541 | } |
---|
1542 | } |
---|
1543 | |
---|
1544 | //printf("Up on %d, status is %d, obj %g its %d cost %g finished %d inf %d infobj %d\n", |
---|
1545 | // choice[i].objectNumber,iStatus,newObjectiveValue,choice[i].numItersUp, |
---|
1546 | // choice[i].upMovement,choice[i].finishedUp,choice[i].numIntInfeasUp, |
---|
1547 | // choice[i].numObjInfeasUp); |
---|
1548 | |
---|
1549 | /* |
---|
1550 | End of evaluation for this candidate variable. Possibilities are: |
---|
1551 | * Both sides below cutoff; this variable is a candidate for branching. |
---|
1552 | * Both sides infeasible or above the objective cutoff: no further action |
---|
1553 | here. Break from the evaluation loop and assume the node will be purged |
---|
1554 | by the caller. |
---|
1555 | * One side below cutoff: Install the branch (i.e., fix the variable). Break |
---|
1556 | from the evaluation loop and assume the node will be reoptimised by the |
---|
1557 | caller. |
---|
1558 | */ |
---|
1559 | if (choice[i].upMovement<1.0e100) { |
---|
1560 | if(choice[i].downMovement<1.0e100) { |
---|
1561 | // feasible - no action |
---|
1562 | } else { |
---|
1563 | // up feasible, down infeasible |
---|
1564 | anyAction=-1; |
---|
1565 | if (!solveAll) { |
---|
1566 | choice[i].possibleBranch->way(1); |
---|
1567 | choice[i].possibleBranch->branch(); |
---|
1568 | break; |
---|
1569 | } else { |
---|
1570 | choice[i].fix=1; |
---|
1571 | } |
---|
1572 | } |
---|
1573 | } else { |
---|
1574 | if(choice[i].downMovement<1.0e100) { |
---|
1575 | // down feasible, up infeasible |
---|
1576 | anyAction=-1; |
---|
1577 | if (!solveAll) { |
---|
1578 | choice[i].possibleBranch->way(-1); |
---|
1579 | choice[i].possibleBranch->branch(); |
---|
1580 | break; |
---|
1581 | } else { |
---|
1582 | choice[i].fix=-1; |
---|
1583 | } |
---|
1584 | } else { |
---|
1585 | // neither side feasible |
---|
1586 | anyAction=-2; |
---|
1587 | break; |
---|
1588 | } |
---|
1589 | } |
---|
1590 | bool hitMaxTime = ( CoinCpuTime()-model->getDblParam(CbcModel::CbcStartSeconds) > |
---|
1591 | model->getDblParam(CbcModel::CbcMaximumSeconds)); |
---|
1592 | if (hitMaxTime) { |
---|
1593 | numberStrong=i+1; |
---|
1594 | break; |
---|
1595 | } |
---|
1596 | } |
---|
1597 | if (!clp) { |
---|
1598 | // Delete the snapshot |
---|
1599 | solver->unmarkHotStart(); |
---|
1600 | } else { |
---|
1601 | delete [] newLower; |
---|
1602 | delete [] newUpper; |
---|
1603 | delete [] outputStuff; |
---|
1604 | int i; |
---|
1605 | for (i=0;i<2*numberStrong;i++) |
---|
1606 | delete [] outputSolution[i]; |
---|
1607 | delete [] outputSolution; |
---|
1608 | } |
---|
1609 | solver->setIntParam(OsiMaxNumIterationHotStart,saveLimit); |
---|
1610 | // restore basis |
---|
1611 | solver->setWarmStart(ws); |
---|
1612 | // Unless infeasible we will carry on |
---|
1613 | // But we could fix anyway |
---|
1614 | if (anyAction==-1&&solveAll) { |
---|
1615 | // apply and take off |
---|
1616 | for (i = 0 ; i < numberStrong ; i++) { |
---|
1617 | if (choice[i].fix) { |
---|
1618 | choice[i].possibleBranch->way(choice[i].fix) ; |
---|
1619 | choice[i].possibleBranch->branch() ; |
---|
1620 | } |
---|
1621 | } |
---|
1622 | bool feasible=true; |
---|
1623 | if (checkFeasibility) { |
---|
1624 | // check branching did not make infeasible |
---|
1625 | int iColumn; |
---|
1626 | int numberColumns = solver->getNumCols(); |
---|
1627 | const double * columnLower = solver->getColLower(); |
---|
1628 | const double * columnUpper = solver->getColUpper(); |
---|
1629 | for (iColumn= 0;iColumn<numberColumns;iColumn++) { |
---|
1630 | if (columnLower[iColumn]>columnUpper[iColumn]+1.0e-5) |
---|
1631 | feasible=false; |
---|
1632 | } |
---|
1633 | } |
---|
1634 | if (feasible) { |
---|
1635 | // can do quick optimality check |
---|
1636 | int easy=2; |
---|
1637 | solver->setHintParam(OsiDoInBranchAndCut,true,OsiHintDo,&easy) ; |
---|
1638 | solver->resolve() ; |
---|
1639 | solver->setHintParam(OsiDoInBranchAndCut,true,OsiHintDo,NULL) ; |
---|
1640 | feasible = solver->isProvenOptimal(); |
---|
1641 | } |
---|
1642 | if (feasible) { |
---|
1643 | memcpy(saveSolution,solver->getColSolution(),numberColumns*sizeof(double)); |
---|
1644 | model->reserveCurrentSolution(saveSolution); |
---|
1645 | memcpy(saveLower,solver->getColLower(),numberColumns*sizeof(double)); |
---|
1646 | memcpy(saveUpper,solver->getColUpper(),numberColumns*sizeof(double)); |
---|
1647 | // Clean up all candidates whih are fixed |
---|
1648 | int numberLeft=0; |
---|
1649 | for (i = 0 ; i < numberStrong ; i++) { |
---|
1650 | Strong thisChoice = choice[i]; |
---|
1651 | choice[i].possibleBranch=NULL; |
---|
1652 | const CbcObject * object = model->object(thisChoice.objectNumber); |
---|
1653 | int preferredWay; |
---|
1654 | double infeasibility = object->infeasibility(preferredWay); |
---|
1655 | if (!infeasibility) { |
---|
1656 | // take out |
---|
1657 | delete thisChoice.possibleBranch; |
---|
1658 | } else { |
---|
1659 | choice[numberLeft++]=thisChoice; |
---|
1660 | } |
---|
1661 | } |
---|
1662 | numberStrong=numberLeft; |
---|
1663 | for (;i<maximumStrong;i++) { |
---|
1664 | delete choice[i].possibleBranch; |
---|
1665 | choice[i].possibleBranch=NULL; |
---|
1666 | } |
---|
1667 | // If all fixed then round again |
---|
1668 | if (!numberLeft) { |
---|
1669 | finished=false; |
---|
1670 | numberStrong=0; |
---|
1671 | saveNumberStrong=0; |
---|
1672 | maximumStrong=1; |
---|
1673 | } else { |
---|
1674 | anyAction=0; |
---|
1675 | } |
---|
1676 | // If these two uncommented then different action |
---|
1677 | anyAction=-1; |
---|
1678 | finished=true; |
---|
1679 | //printf("some fixed but continuing %d left\n",numberLeft); |
---|
1680 | } else { |
---|
1681 | anyAction=-2; // say infeasible |
---|
1682 | } |
---|
1683 | } |
---|
1684 | delete ws; |
---|
1685 | |
---|
1686 | /* |
---|
1687 | anyAction >= 0 indicates that strong branching didn't produce any monotone |
---|
1688 | variables. Sift through the candidates for the best one. |
---|
1689 | |
---|
1690 | QUERY: Setting numberNodes looks to be a distributed noop. numberNodes is |
---|
1691 | local to this code block. Perhaps should be numberNodes_ from model? |
---|
1692 | Unclear what this calculation is doing. |
---|
1693 | */ |
---|
1694 | if (anyAction>=0) { |
---|
1695 | |
---|
1696 | int numberNodes = model->getNodeCount(); |
---|
1697 | // get average cost per iteration and assume stopped ones |
---|
1698 | // would stop after 50% more iterations at average cost??? !!! ??? |
---|
1699 | double averageCostPerIteration=0.0; |
---|
1700 | double totalNumberIterations=1.0; |
---|
1701 | int smallestNumberInfeasibilities=INT_MAX; |
---|
1702 | for (i=0;i<numberStrong;i++) { |
---|
1703 | totalNumberIterations += choice[i].numItersDown + |
---|
1704 | choice[i].numItersUp ; |
---|
1705 | averageCostPerIteration += choice[i].downMovement + |
---|
1706 | choice[i].upMovement; |
---|
1707 | smallestNumberInfeasibilities= |
---|
1708 | CoinMin(CoinMin(choice[i].numIntInfeasDown , |
---|
1709 | choice[i].numIntInfeasUp ), |
---|
1710 | smallestNumberInfeasibilities); |
---|
1711 | } |
---|
1712 | if (smallestNumberInfeasibilities>=numberIntegerInfeasibilities) |
---|
1713 | numberNodes=1000000; // switch off search for better solution |
---|
1714 | numberNodes=1000000; // switch off anyway |
---|
1715 | averageCostPerIteration /= totalNumberIterations; |
---|
1716 | // all feasible - choose best bet |
---|
1717 | |
---|
1718 | #if 0 |
---|
1719 | for (i = 0 ; i < numberStrong ; i++) |
---|
1720 | { int iColumn = |
---|
1721 | model->integerVariable()[choice[i].possibleBranch->variable()] ; |
---|
1722 | model->messageHandler()->message(CBC_STRONG,model->messages()) |
---|
1723 | << i << iColumn |
---|
1724 | <<choice[i].downMovement<<choice[i].numIntInfeasDown |
---|
1725 | <<choice[i].upMovement<<choice[i].numIntInfeasUp |
---|
1726 | <<choice[i].possibleBranch->value() |
---|
1727 | <<CoinMessageEol; |
---|
1728 | int betterWay = decision->betterBranch(choice[i].possibleBranch, |
---|
1729 | branch_, |
---|
1730 | choice[i].upMovement, |
---|
1731 | choice[i].numIntInfeasUp , |
---|
1732 | choice[i].downMovement, |
---|
1733 | choice[i].numIntInfeasDown ); |
---|
1734 | if (betterWay) { |
---|
1735 | delete branch_; |
---|
1736 | // C) create branching object |
---|
1737 | branch_ = choice[i].possibleBranch; |
---|
1738 | choice[i].possibleBranch=NULL; |
---|
1739 | branch_->way(betterWay); |
---|
1740 | } |
---|
1741 | } |
---|
1742 | #else |
---|
1743 | // New method does all at once so it can be more sophisticated |
---|
1744 | // in deciding how to balance actions. |
---|
1745 | // But it does need arrays |
---|
1746 | double * changeUp = new double [numberStrong]; |
---|
1747 | int * numberInfeasibilitiesUp = new int [numberStrong]; |
---|
1748 | double * changeDown = new double [numberStrong]; |
---|
1749 | int * numberInfeasibilitiesDown = new int [numberStrong]; |
---|
1750 | CbcBranchingObject ** objects = new CbcBranchingObject * [ numberStrong]; |
---|
1751 | for (i = 0 ; i < numberStrong ; i++) { |
---|
1752 | int iColumn = choice[i].possibleBranch->variable() ; |
---|
1753 | model->messageHandler()->message(CBC_STRONG,model->messages()) |
---|
1754 | << i << iColumn |
---|
1755 | <<choice[i].downMovement<<choice[i].numIntInfeasDown |
---|
1756 | <<choice[i].upMovement<<choice[i].numIntInfeasUp |
---|
1757 | <<choice[i].possibleBranch->value() |
---|
1758 | <<CoinMessageEol; |
---|
1759 | changeUp[i]=choice[i].upMovement; |
---|
1760 | numberInfeasibilitiesUp[i] = choice[i].numIntInfeasUp; |
---|
1761 | changeDown[i]=choice[i].downMovement; |
---|
1762 | numberInfeasibilitiesDown[i] = choice[i].numIntInfeasDown; |
---|
1763 | objects[i] = choice[i].possibleBranch; |
---|
1764 | } |
---|
1765 | int whichObject = decision->bestBranch(objects,numberStrong,numberUnsatisfied_, |
---|
1766 | changeUp,numberInfeasibilitiesUp, |
---|
1767 | changeDown,numberInfeasibilitiesDown, |
---|
1768 | objectiveValue); |
---|
1769 | // move branching object and make sure it will not be deleted |
---|
1770 | if (whichObject>=0) { |
---|
1771 | branch_ = objects[whichObject]; |
---|
1772 | choice[whichObject].possibleBranch=NULL; |
---|
1773 | } |
---|
1774 | delete [] changeUp; |
---|
1775 | delete [] numberInfeasibilitiesUp; |
---|
1776 | delete [] changeDown; |
---|
1777 | delete [] numberInfeasibilitiesDown; |
---|
1778 | delete [] objects; |
---|
1779 | #endif |
---|
1780 | } |
---|
1781 | } |
---|
1782 | /* |
---|
1783 | Simple branching. Probably just one, but we may have got here |
---|
1784 | because of an odd branch e.g. a cut |
---|
1785 | */ |
---|
1786 | else { |
---|
1787 | // not strong |
---|
1788 | // C) create branching object |
---|
1789 | branch_ = choice[bestChoice].possibleBranch; |
---|
1790 | choice[bestChoice].possibleBranch=NULL; |
---|
1791 | } |
---|
1792 | } |
---|
1793 | // Set guessed solution value |
---|
1794 | guessedObjectiveValue_ = objectiveValue_+estimatedDegradation; |
---|
1795 | /* |
---|
1796 | Cleanup, then we're outta here. |
---|
1797 | */ |
---|
1798 | if (!model->branchingMethod()) |
---|
1799 | delete decision; |
---|
1800 | |
---|
1801 | for (i=0;i<maximumStrong;i++) |
---|
1802 | delete choice[i].possibleBranch; |
---|
1803 | delete [] choice; |
---|
1804 | delete [] saveLower; |
---|
1805 | delete [] saveUpper; |
---|
1806 | |
---|
1807 | // restore solution |
---|
1808 | solver->setColSolution(saveSolution); |
---|
1809 | delete [] saveSolution; |
---|
1810 | return anyAction; |
---|
1811 | } |
---|
1812 | |
---|
1813 | |
---|
1814 | CbcNode::CbcNode(const CbcNode & rhs) |
---|
1815 | { |
---|
1816 | #ifdef CHECK_NODE |
---|
1817 | printf("CbcNode %x Constructor from rhs %x\n",this,&rhs); |
---|
1818 | #endif |
---|
1819 | if (rhs.nodeInfo_) |
---|
1820 | nodeInfo_ = rhs.nodeInfo_->clone(); |
---|
1821 | else |
---|
1822 | nodeInfo_=NULL; |
---|
1823 | objectiveValue_=rhs.objectiveValue_; |
---|
1824 | guessedObjectiveValue_ = rhs.guessedObjectiveValue_; |
---|
1825 | if (rhs.branch_) |
---|
1826 | branch_=rhs.branch_->clone(); |
---|
1827 | else |
---|
1828 | branch_=NULL; |
---|
1829 | depth_ = rhs.depth_; |
---|
1830 | numberUnsatisfied_ = rhs.numberUnsatisfied_; |
---|
1831 | nodeNumber_=rhs.nodeNumber_; |
---|
1832 | } |
---|
1833 | |
---|
1834 | CbcNode & |
---|
1835 | CbcNode::operator=(const CbcNode & rhs) |
---|
1836 | { |
---|
1837 | if (this != &rhs) { |
---|
1838 | delete nodeInfo_; |
---|
1839 | if (nodeInfo_) |
---|
1840 | nodeInfo_ = rhs.nodeInfo_->clone(); |
---|
1841 | else |
---|
1842 | nodeInfo_ = NULL; |
---|
1843 | objectiveValue_=rhs.objectiveValue_; |
---|
1844 | guessedObjectiveValue_ = rhs.guessedObjectiveValue_; |
---|
1845 | if (rhs.branch_) |
---|
1846 | branch_=rhs.branch_->clone(); |
---|
1847 | else |
---|
1848 | branch_=NULL, |
---|
1849 | depth_ = rhs.depth_; |
---|
1850 | numberUnsatisfied_ = rhs.numberUnsatisfied_; |
---|
1851 | nodeNumber_=rhs.nodeNumber_; |
---|
1852 | } |
---|
1853 | return *this; |
---|
1854 | } |
---|
1855 | |
---|
1856 | |
---|
1857 | CbcNode::~CbcNode () |
---|
1858 | { |
---|
1859 | #ifdef CHECK_NODE |
---|
1860 | if (nodeInfo_) |
---|
1861 | printf("CbcNode %x Destructor nodeInfo %x (%d)\n", |
---|
1862 | this,nodeInfo_,nodeInfo_->numberPointingToThis()); |
---|
1863 | else |
---|
1864 | printf("CbcNode %x Destructor nodeInfo %x (?)\n", |
---|
1865 | this,nodeInfo_); |
---|
1866 | #endif |
---|
1867 | if (nodeInfo_) { |
---|
1868 | nodeInfo_->nullOwner(); |
---|
1869 | int numberToDelete=nodeInfo_->numberBranchesLeft(); |
---|
1870 | // CbcNodeInfo * parent = nodeInfo_->parent(); |
---|
1871 | if (nodeInfo_->decrement(numberToDelete)==0) { |
---|
1872 | delete nodeInfo_; |
---|
1873 | } else { |
---|
1874 | //printf("node %x nodeinfo %x parent %x\n",this,nodeInfo_,parent); |
---|
1875 | // anyway decrement parent |
---|
1876 | //if (parent) |
---|
1877 | ///parent->decrement(1); |
---|
1878 | } |
---|
1879 | } |
---|
1880 | delete branch_; |
---|
1881 | } |
---|
1882 | // Decrement active cut counts |
---|
1883 | void |
---|
1884 | CbcNode::decrementCuts(int change) |
---|
1885 | { |
---|
1886 | if(nodeInfo_) { |
---|
1887 | nodeInfo_->decrementCuts(change); |
---|
1888 | } |
---|
1889 | } |
---|
1890 | void |
---|
1891 | CbcNode::decrementParentCuts(int change) |
---|
1892 | { |
---|
1893 | if(nodeInfo_) { |
---|
1894 | nodeInfo_->decrementParentCuts(change); |
---|
1895 | } |
---|
1896 | } |
---|
1897 | |
---|
1898 | /* |
---|
1899 | Initialize reference counts (numberPointingToThis, numberBranchesLeft_) |
---|
1900 | in the attached nodeInfo_. |
---|
1901 | */ |
---|
1902 | void |
---|
1903 | CbcNode::initializeInfo() |
---|
1904 | { |
---|
1905 | assert(nodeInfo_ && branch_) ; |
---|
1906 | nodeInfo_->initializeInfo(branch_->numberBranches()); |
---|
1907 | } |
---|
1908 | // Nulls out node info |
---|
1909 | void |
---|
1910 | CbcNode::nullNodeInfo() |
---|
1911 | { |
---|
1912 | nodeInfo_=NULL; |
---|
1913 | } |
---|
1914 | |
---|
1915 | int |
---|
1916 | CbcNode::branch() |
---|
1917 | { |
---|
1918 | double changeInGuessed=branch_->branch(true); |
---|
1919 | guessedObjectiveValue_+= changeInGuessed; |
---|
1920 | return nodeInfo_->branchedOn(); |
---|
1921 | } |
---|