1 | /********************************************************************** |
---|
2 | TODO: |
---|
3 | |
---|
4 | If we swap two nodes (parent and grandparent) then we should check if anything |
---|
5 | has been explored under GP after the swap, and if not, just get rid of GP and |
---|
6 | everything below. |
---|
7 | |
---|
8 | If strong branching fixed anything in the grandparent we may still swap it |
---|
9 | with the parent (we don't do it yet, see the test on |
---|
10 | strong_branching_fixed_vars_), but those fixings must be moved to the parent as |
---|
11 | well. |
---|
12 | |
---|
13 | Same thing for locally valid cuts, if GP has them and GP is switched with P |
---|
14 | then locally valid cuts must be treated as generated in P. |
---|
15 | |
---|
16 | Same for reduced cost fixing :-(. If P has done any then P and GP can't be |
---|
17 | switched. |
---|
18 | |
---|
19 | Maybe the best solution is to create local information (strong branching, rc |
---|
20 | fixing, local cuts, etc.) only every so often, say every 10th level. Then that |
---|
21 | would block switches, but everywhere else we would be safe. Good question |
---|
22 | what is worth more: the switches or the local info. |
---|
23 | |
---|
24 | Alternative solution: Do not add local info to the tree, but keep it in a set |
---|
25 | of excluded clauses ala Todias Achterberg: Conflict Analysis in Mixed Integer |
---|
26 | Programming. |
---|
27 | |
---|
28 | We may want to disable fixing by strong branching completely and if such |
---|
29 | situation occurs then simply do the branch and one side will be fathomed |
---|
30 | immediately and we can try to switch. |
---|
31 | |
---|
32 | Bound modifications when nodes are swapped could be avoided if always start |
---|
33 | from original bounds and go back to root to apply all branching decisions. On |
---|
34 | the other hand, reduced cost fixing would be lost. And so would fixings by |
---|
35 | strong branching. Although both could be stored in an array keeping track of |
---|
36 | changes implied by the branching decisions. |
---|
37 | |
---|
38 | ********************************************************************** |
---|
39 | |
---|
40 | |
---|
41 | **********************************************************************/ |
---|
42 | #include "CoinTime.hpp" |
---|
43 | #include "OsiClpSolverInterface.hpp" |
---|
44 | |
---|
45 | #define FUNNY_BRANCHING 1 |
---|
46 | #define DEBUG_DYNAMIC_BRANCHING |
---|
47 | |
---|
48 | #ifdef DEBUG_DYNAMIC_BRANCHING |
---|
49 | int dyn_debug = 1; |
---|
50 | #endif |
---|
51 | |
---|
52 | // below needed for pathetic branch and bound code |
---|
53 | #include <vector> |
---|
54 | #include <map> |
---|
55 | |
---|
56 | using namespace std; |
---|
57 | |
---|
58 | class DBVectorNode; |
---|
59 | |
---|
60 | class LPresult { |
---|
61 | private: |
---|
62 | LPresult(const LPresult& rhs); |
---|
63 | LPresult& operator=(const LPresult& rhs); |
---|
64 | void gutsOfConstructor(const OsiSolverInterface& model); |
---|
65 | public: |
---|
66 | LPresult(const OsiSolverInterface& model); |
---|
67 | ~LPresult(); |
---|
68 | public: |
---|
69 | bool isAbandoned; |
---|
70 | bool isProvenDualInfeasible; |
---|
71 | bool isPrimalObjectiveLimitReached; |
---|
72 | bool isProvenOptimal; |
---|
73 | bool isDualObjectiveLimitReached; |
---|
74 | bool isIterationLimitReached; |
---|
75 | bool isProvenPrimalInfeasible; |
---|
76 | double getObjSense; |
---|
77 | double getObjValue; |
---|
78 | double* getReducedCost; |
---|
79 | double* getColLower; |
---|
80 | double* getColUpper; |
---|
81 | double* rowDualRay; |
---|
82 | double* colDualRay; |
---|
83 | double* getColSolution; // FIXME Not needed, just for debugging |
---|
84 | double yb_plus_rl_minus_su; |
---|
85 | }; |
---|
86 | |
---|
87 | void |
---|
88 | LPresult::gutsOfConstructor(const OsiSolverInterface& model) |
---|
89 | { |
---|
90 | isAbandoned = model.isAbandoned(); |
---|
91 | isProvenDualInfeasible = model.isProvenDualInfeasible(); |
---|
92 | isPrimalObjectiveLimitReached = model.isPrimalObjectiveLimitReached(); |
---|
93 | isProvenOptimal = model.isProvenOptimal(); |
---|
94 | isDualObjectiveLimitReached = model.isDualObjectiveLimitReached(); |
---|
95 | isIterationLimitReached = model.isIterationLimitReached(); |
---|
96 | isProvenPrimalInfeasible = model.isProvenPrimalInfeasible(); |
---|
97 | getObjSense = model.getObjSense(); |
---|
98 | getObjValue = model.getObjValue(); |
---|
99 | |
---|
100 | getReducedCost = new double[model.getNumCols()]; |
---|
101 | CoinDisjointCopyN(model.getReducedCost(), model.getNumCols(), getReducedCost); |
---|
102 | getColLower = new double[model.getNumCols()]; |
---|
103 | CoinDisjointCopyN(model.getColLower(), model.getNumCols(), getColLower); |
---|
104 | getColUpper = new double[model.getNumCols()]; |
---|
105 | CoinDisjointCopyN(model.getColUpper(), model.getNumCols(), getColUpper); |
---|
106 | getColSolution = new double[model.getNumCols()]; |
---|
107 | CoinDisjointCopyN(model.getColSolution(), model.getNumCols(), getColSolution); |
---|
108 | rowDualRay = NULL; |
---|
109 | colDualRay = NULL; |
---|
110 | |
---|
111 | yb_plus_rl_minus_su = 0; |
---|
112 | |
---|
113 | if (isIterationLimitReached || |
---|
114 | isAbandoned || |
---|
115 | isProvenDualInfeasible || |
---|
116 | isPrimalObjectiveLimitReached) { |
---|
117 | return; // None of these should happen |
---|
118 | } |
---|
119 | // If isDualObjectiveLimitReached is true then isProvenOptimal and |
---|
120 | // isProvenPrimalInfeasible are correct (search for mustResolve). |
---|
121 | // If isDualObjectiveLimitReached false, then the others were correct even |
---|
122 | // without resolve. So... we need some data only if we are infeasible, but |
---|
123 | // that flag is correct, so we can test it. |
---|
124 | if (isProvenPrimalInfeasible) { |
---|
125 | const std::vector<double*> dual_rays = model.getDualRays(1); |
---|
126 | if (dual_rays.size() == 0) { |
---|
127 | printf("WARNING: LP is infeas, but no dual ray is returned!\n"); |
---|
128 | return; |
---|
129 | } |
---|
130 | const double* rub = model.getRowUpper(); |
---|
131 | const double* rlb = model.getRowLower(); |
---|
132 | const double* cub = model.getColUpper(); |
---|
133 | const double* clb = model.getColLower(); |
---|
134 | const int numRows = model.getNumRows(); |
---|
135 | const int numCols = model.getNumCols(); |
---|
136 | rowDualRay = dual_rays[0]; |
---|
137 | colDualRay = new double[model.getNumCols()]; |
---|
138 | #if 1 |
---|
139 | // WARNING ==== WARNING ==== WARNING ==== WARNING ==== WARNING |
---|
140 | // For some reason clp returns the negative of the dual ray... |
---|
141 | for (int i = 0; i < numRows; ++i) { |
---|
142 | rowDualRay[i] = -rowDualRay[i]; |
---|
143 | } |
---|
144 | #endif |
---|
145 | const CoinPackedMatrix* m = model.getMatrixByCol(); |
---|
146 | m->transposeTimes(rowDualRay, colDualRay); |
---|
147 | |
---|
148 | for (int i = 0; i < numRows; ++i) { |
---|
149 | const double ray_i = rowDualRay[i]; |
---|
150 | if (ray_i > 1e-8) { |
---|
151 | yb_plus_rl_minus_su += ray_i*rlb[i]; |
---|
152 | assert(rlb[i] > -1e29); |
---|
153 | } else if (ray_i < 1e-8) { |
---|
154 | yb_plus_rl_minus_su += ray_i*rub[i]; |
---|
155 | assert(rub[i] < 1e29); |
---|
156 | } |
---|
157 | } |
---|
158 | for (int j = 0; j < numCols; ++j) { |
---|
159 | if (colDualRay[j] > 1e-8) { // |
---|
160 | yb_plus_rl_minus_su -= colDualRay[j] * cub[j]; |
---|
161 | assert(cub[j] < 1e29); |
---|
162 | } else if (colDualRay[j] < -1e-8) { |
---|
163 | yb_plus_rl_minus_su -= colDualRay[j] * clb[j]; |
---|
164 | assert(clb[j] > -1e29); |
---|
165 | } |
---|
166 | } |
---|
167 | for (int i = dual_rays.size()-1; i > 0; --i) { |
---|
168 | delete[] dual_rays[i]; |
---|
169 | } |
---|
170 | assert(yb_plus_rl_minus_su > 1e-8); |
---|
171 | } |
---|
172 | } |
---|
173 | |
---|
174 | LPresult::LPresult(const OsiSolverInterface& model) |
---|
175 | { |
---|
176 | gutsOfConstructor(model); |
---|
177 | } |
---|
178 | |
---|
179 | LPresult::~LPresult() |
---|
180 | { |
---|
181 | delete[] getReducedCost; |
---|
182 | delete[] getColLower; |
---|
183 | delete[] getColUpper; |
---|
184 | delete[] getColSolution; |
---|
185 | delete[] colDualRay; |
---|
186 | delete[] rowDualRay; |
---|
187 | } |
---|
188 | |
---|
189 | // Trivial class for Branch and Bound |
---|
190 | |
---|
191 | class DBNodeSimple { |
---|
192 | public: |
---|
193 | enum DBNodeWay { |
---|
194 | WAY_UNSET=0x00, |
---|
195 | WAY_DOWN_UP__NOTHING_DONE=0x10, |
---|
196 | WAY_DOWN_UP__DOWN_DONE=0x11, |
---|
197 | WAY_DOWN_UP__BOTH_DONE=0x13, |
---|
198 | WAY_DOWN_UP__=0x10, |
---|
199 | WAY_UP_DOWN__NOTHING_DONE=0x20, |
---|
200 | WAY_UP_DOWN__UP_DONE=0x22, |
---|
201 | WAY_UP_DOWN__BOTH_DONE=0x23, |
---|
202 | WAY_UP_DOWN__=0x20, |
---|
203 | WAY_DOWN_CURRENT=0x01, |
---|
204 | WAY_UP_CURRENT=0x02, |
---|
205 | WAY_BOTH_DONE=0x03, |
---|
206 | |
---|
207 | }; |
---|
208 | |
---|
209 | private: |
---|
210 | void gutsOfCopy(const DBNodeSimple& rhs); |
---|
211 | void gutsOfConstructor (int index, OsiSolverInterface &model, |
---|
212 | int numberIntegers, int * integer, |
---|
213 | CoinWarmStart * basis); |
---|
214 | void gutsOfDestructor(); |
---|
215 | |
---|
216 | public: |
---|
217 | |
---|
218 | // Default Constructor |
---|
219 | DBNodeSimple (); |
---|
220 | |
---|
221 | // Constructor from current state (and list of integers) |
---|
222 | // Also chooses branching variable (if none set to -1) |
---|
223 | DBNodeSimple (int index, OsiSolverInterface &model, |
---|
224 | int numberIntegers, int * integer, |
---|
225 | CoinWarmStart * basis); |
---|
226 | // Copy constructor |
---|
227 | DBNodeSimple ( const DBNodeSimple &); |
---|
228 | |
---|
229 | // Assignment operator |
---|
230 | DBNodeSimple & operator=( const DBNodeSimple& rhs); |
---|
231 | |
---|
232 | // Destructor |
---|
233 | ~DBNodeSimple (); |
---|
234 | |
---|
235 | // Extension - true if other extension of this |
---|
236 | bool extension(const DBNodeSimple & other, |
---|
237 | const double * originalLower, |
---|
238 | const double * originalUpper) const; |
---|
239 | // Tests if we can switch this node (this is the parent) with its parent |
---|
240 | bool canSwitchParentWithGrandparent(const int* which, |
---|
241 | LPresult& lpres, |
---|
242 | const int * original_lower, |
---|
243 | const int * original_upper, |
---|
244 | DBVectorNode & branchingTree); |
---|
245 | inline bool bothChildDone() const { |
---|
246 | return (way_ & WAY_BOTH_DONE) == WAY_BOTH_DONE; |
---|
247 | } |
---|
248 | inline bool workingOnDownChild() const { |
---|
249 | return (way_ == WAY_DOWN_UP__DOWN_DONE || way_ == WAY_UP_DOWN__BOTH_DONE); |
---|
250 | } |
---|
251 | /* Return whether the child that will be now processed is the down branch or |
---|
252 | not. */ |
---|
253 | inline bool advanceWay() { |
---|
254 | switch (way_) { |
---|
255 | case WAY_UNSET: |
---|
256 | case WAY_DOWN_UP__BOTH_DONE: |
---|
257 | case WAY_UP_DOWN__BOTH_DONE: |
---|
258 | abort(); |
---|
259 | case WAY_DOWN_UP__NOTHING_DONE: |
---|
260 | way_ = WAY_DOWN_UP__DOWN_DONE; |
---|
261 | return true; |
---|
262 | case WAY_DOWN_UP__DOWN_DONE: |
---|
263 | way_ = WAY_DOWN_UP__BOTH_DONE; |
---|
264 | return false; |
---|
265 | case WAY_UP_DOWN__NOTHING_DONE: |
---|
266 | way_ = WAY_UP_DOWN__UP_DONE; |
---|
267 | return false; |
---|
268 | case WAY_UP_DOWN__UP_DONE: |
---|
269 | way_ = WAY_UP_DOWN__BOTH_DONE; |
---|
270 | return true; |
---|
271 | } |
---|
272 | return true; // to placate some compilers... |
---|
273 | } |
---|
274 | |
---|
275 | |
---|
276 | // Public data |
---|
277 | int index_; |
---|
278 | // The id of the node |
---|
279 | int node_id_; |
---|
280 | // Basis (should use tree, but not as wasteful as bounds!) |
---|
281 | CoinWarmStart * basis_; |
---|
282 | // Objective value (COIN_DBL_MAX) if spare node |
---|
283 | double objectiveValue_; |
---|
284 | // Branching variable (0 is first integer) |
---|
285 | int variable_; |
---|
286 | // Way to branch: see enum DBNodeWay |
---|
287 | int way_; |
---|
288 | // Number of integers (for length of arrays) |
---|
289 | int numberIntegers_; |
---|
290 | // Current value |
---|
291 | double value_; |
---|
292 | // Parent |
---|
293 | int parent_; |
---|
294 | // Left child |
---|
295 | int child_down_; |
---|
296 | // Right child |
---|
297 | int child_up_; |
---|
298 | // Whether strong branching fixed variables when we branched on this node |
---|
299 | bool strong_branching_fixed_vars_; |
---|
300 | // Whether reduced cost fixing fixed anything in this node |
---|
301 | bool reduced_cost_fixed_vars_; |
---|
302 | // Previous in chain |
---|
303 | int previous_; |
---|
304 | // Next in chain |
---|
305 | int next_; |
---|
306 | // Now I must use tree |
---|
307 | // Bounds stored in full (for integers) |
---|
308 | int * lower_; |
---|
309 | int * upper_; |
---|
310 | }; |
---|
311 | |
---|
312 | |
---|
313 | DBNodeSimple::DBNodeSimple() : |
---|
314 | index_(-1), |
---|
315 | node_id_(-1), |
---|
316 | basis_(NULL), |
---|
317 | objectiveValue_(COIN_DBL_MAX), |
---|
318 | variable_(-100), |
---|
319 | way_(WAY_UNSET), |
---|
320 | numberIntegers_(0), |
---|
321 | value_(0.5), |
---|
322 | parent_(-1), |
---|
323 | child_down_(-1), |
---|
324 | child_up_(-1), |
---|
325 | strong_branching_fixed_vars_(false), |
---|
326 | reduced_cost_fixed_vars_(false), |
---|
327 | previous_(-1), |
---|
328 | next_(-1), |
---|
329 | lower_(NULL), |
---|
330 | upper_(NULL) |
---|
331 | { |
---|
332 | } |
---|
333 | DBNodeSimple::DBNodeSimple(int index, OsiSolverInterface & model, |
---|
334 | int numberIntegers, int * integer,CoinWarmStart * basis) |
---|
335 | { |
---|
336 | gutsOfConstructor(index, model,numberIntegers,integer,basis); |
---|
337 | } |
---|
338 | void |
---|
339 | DBNodeSimple::gutsOfConstructor(int index, OsiSolverInterface & model, |
---|
340 | int numberIntegers, int * integer,CoinWarmStart * basis) |
---|
341 | { |
---|
342 | index_ = index; |
---|
343 | node_id_ = -1; |
---|
344 | basis_ = basis; |
---|
345 | variable_=-1; |
---|
346 | way_ = WAY_UNSET; |
---|
347 | numberIntegers_=numberIntegers; |
---|
348 | value_=0.0; |
---|
349 | parent_ = -1; |
---|
350 | child_down_ = -1; |
---|
351 | child_up_ = -1; |
---|
352 | strong_branching_fixed_vars_ = false; |
---|
353 | reduced_cost_fixed_vars_ = false; |
---|
354 | previous_ = -1; |
---|
355 | next_ = -1; |
---|
356 | if (model.isProvenOptimal()&&!model.isDualObjectiveLimitReached()) { |
---|
357 | objectiveValue_ = model.getObjSense()*model.getObjValue(); |
---|
358 | } else { |
---|
359 | objectiveValue_ = 1.0e100; |
---|
360 | lower_ = NULL; |
---|
361 | upper_ = NULL; |
---|
362 | return; // node cutoff |
---|
363 | } |
---|
364 | lower_ = new int [numberIntegers_]; |
---|
365 | upper_ = new int [numberIntegers_]; |
---|
366 | assert (upper_!=NULL); |
---|
367 | const double * lower = model.getColLower(); |
---|
368 | const double * upper = model.getColUpper(); |
---|
369 | const double * solution = model.getColSolution(); |
---|
370 | int i; |
---|
371 | // Hard coded integer tolerance |
---|
372 | #define INTEGER_TOLERANCE 1.0e-6 |
---|
373 | ///////// Start of Strong branching code - can be ignored |
---|
374 | // Number of strong branching candidates |
---|
375 | #define STRONG_BRANCHING 5 |
---|
376 | #ifdef STRONG_BRANCHING |
---|
377 | double upMovement[STRONG_BRANCHING]; |
---|
378 | double downMovement[STRONG_BRANCHING]; |
---|
379 | double solutionValue[STRONG_BRANCHING]; |
---|
380 | int chosen[STRONG_BRANCHING]; |
---|
381 | int iSmallest=0; |
---|
382 | // initialize distance from integer |
---|
383 | for (i=0;i<STRONG_BRANCHING;i++) { |
---|
384 | upMovement[i]=0.0; |
---|
385 | chosen[i]=-1; |
---|
386 | } |
---|
387 | variable_=-1; |
---|
388 | // This has hard coded integer tolerance |
---|
389 | double mostAway=INTEGER_TOLERANCE; |
---|
390 | int numberAway=0; |
---|
391 | for (i=0;i<numberIntegers;i++) { |
---|
392 | int iColumn = integer[i]; |
---|
393 | lower_[i]=(int)lower[iColumn]; |
---|
394 | upper_[i]=(int)upper[iColumn]; |
---|
395 | double value = solution[iColumn]; |
---|
396 | value = max(value,(double) lower_[i]); |
---|
397 | value = min(value,(double) upper_[i]); |
---|
398 | double nearest = floor(value+0.5); |
---|
399 | if (fabs(value-nearest)>INTEGER_TOLERANCE) |
---|
400 | numberAway++; |
---|
401 | if (fabs(value-nearest)>mostAway) { |
---|
402 | double away = fabs(value-nearest); |
---|
403 | if (away>upMovement[iSmallest]) { |
---|
404 | //add to list |
---|
405 | upMovement[iSmallest]=away; |
---|
406 | solutionValue[iSmallest]=value; |
---|
407 | chosen[iSmallest]=i; |
---|
408 | int j; |
---|
409 | iSmallest=-1; |
---|
410 | double smallest = 1.0; |
---|
411 | for (j=0;j<STRONG_BRANCHING;j++) { |
---|
412 | if (upMovement[j]<smallest) { |
---|
413 | smallest=upMovement[j]; |
---|
414 | iSmallest=j; |
---|
415 | } |
---|
416 | } |
---|
417 | } |
---|
418 | } |
---|
419 | } |
---|
420 | int numberStrong=0; |
---|
421 | for (i=0;i<STRONG_BRANCHING;i++) { |
---|
422 | if (chosen[i]>=0) { |
---|
423 | numberStrong ++; |
---|
424 | variable_ = chosen[i]; |
---|
425 | } |
---|
426 | } |
---|
427 | // out strong branching if bit set |
---|
428 | OsiClpSolverInterface* clp = |
---|
429 | dynamic_cast<OsiClpSolverInterface*>(&model); |
---|
430 | if (clp&&(clp->specialOptions()&16)!=0&&numberStrong>1) { |
---|
431 | int j; |
---|
432 | int iBest=-1; |
---|
433 | double best = 0.0; |
---|
434 | for (j=0;j<STRONG_BRANCHING;j++) { |
---|
435 | if (upMovement[j]>best) { |
---|
436 | best=upMovement[j]; |
---|
437 | iBest=j; |
---|
438 | } |
---|
439 | } |
---|
440 | numberStrong=1; |
---|
441 | variable_=chosen[iBest]; |
---|
442 | } |
---|
443 | if (numberStrong==1) { |
---|
444 | // just one - makes it easy |
---|
445 | int iColumn = integer[variable_]; |
---|
446 | double value = solution[iColumn]; |
---|
447 | value = max(value,(double) lower_[variable_]); |
---|
448 | value = min(value,(double) upper_[variable_]); |
---|
449 | double nearest = floor(value+0.5); |
---|
450 | value_=value; |
---|
451 | if (value<=nearest) |
---|
452 | way_ = WAY_UP_DOWN__NOTHING_DONE; // up |
---|
453 | else |
---|
454 | way_ = WAY_DOWN_UP__NOTHING_DONE; // down |
---|
455 | } else if (numberStrong) { |
---|
456 | // more than one - choose |
---|
457 | bool chooseOne=true; |
---|
458 | model.markHotStart(); |
---|
459 | for (i=0;i<STRONG_BRANCHING;i++) { |
---|
460 | int iInt = chosen[i]; |
---|
461 | if (iInt>=0) { |
---|
462 | int iColumn = integer[iInt]; |
---|
463 | double value = solutionValue[i]; // value of variable in original |
---|
464 | double objectiveChange; |
---|
465 | value = max(value,(double) lower_[iInt]); |
---|
466 | value = min(value,(double) upper_[iInt]); |
---|
467 | |
---|
468 | // try down |
---|
469 | |
---|
470 | model.setColUpper(iColumn,floor(value)); |
---|
471 | model.solveFromHotStart(); |
---|
472 | model.setColUpper(iColumn,upper_[iInt]); |
---|
473 | if (model.isProvenOptimal()&&!model.isDualObjectiveLimitReached()) { |
---|
474 | objectiveChange = model.getObjSense()*model.getObjValue() |
---|
475 | - objectiveValue_; |
---|
476 | } else { |
---|
477 | objectiveChange = 1.0e100; |
---|
478 | } |
---|
479 | assert (objectiveChange>-1.0e-5); |
---|
480 | objectiveChange = CoinMax(objectiveChange,0.0); |
---|
481 | downMovement[i]=objectiveChange; |
---|
482 | |
---|
483 | // try up |
---|
484 | |
---|
485 | model.setColLower(iColumn,ceil(value)); |
---|
486 | model.solveFromHotStart(); |
---|
487 | model.setColLower(iColumn,lower_[iInt]); |
---|
488 | if (model.isProvenOptimal()&&!model.isDualObjectiveLimitReached()) { |
---|
489 | objectiveChange = model.getObjSense()*model.getObjValue() |
---|
490 | - objectiveValue_; |
---|
491 | } else { |
---|
492 | objectiveChange = 1.0e100; |
---|
493 | } |
---|
494 | assert (objectiveChange>-1.0e-5); |
---|
495 | objectiveChange = CoinMax(objectiveChange,0.0); |
---|
496 | upMovement[i]=objectiveChange; |
---|
497 | |
---|
498 | /* Possibilities are: |
---|
499 | Both sides feasible - store |
---|
500 | Neither side feasible - set objective high and exit |
---|
501 | One side feasible - change bounds and resolve |
---|
502 | */ |
---|
503 | bool solveAgain=false; |
---|
504 | if (upMovement[i]<1.0e100) { |
---|
505 | if(downMovement[i]<1.0e100) { |
---|
506 | // feasible - no action |
---|
507 | } else { |
---|
508 | // up feasible, down infeasible |
---|
509 | solveAgain = true; |
---|
510 | model.setColLower(iColumn,ceil(value)); |
---|
511 | } |
---|
512 | } else { |
---|
513 | if(downMovement[i]<1.0e100) { |
---|
514 | // down feasible, up infeasible |
---|
515 | solveAgain = true; |
---|
516 | model.setColUpper(iColumn,floor(value)); |
---|
517 | } else { |
---|
518 | // neither side feasible |
---|
519 | objectiveValue_=1.0e100; |
---|
520 | chooseOne=false; |
---|
521 | break; |
---|
522 | } |
---|
523 | } |
---|
524 | if (solveAgain) { |
---|
525 | // need to solve problem again - signal this |
---|
526 | variable_ = numberIntegers; |
---|
527 | chooseOne=false; |
---|
528 | break; |
---|
529 | } |
---|
530 | } |
---|
531 | } |
---|
532 | if (chooseOne) { |
---|
533 | // choose the one that makes most difference both ways |
---|
534 | double best = -1.0; |
---|
535 | double best2 = -1.0; |
---|
536 | for (i=0;i<STRONG_BRANCHING;i++) { |
---|
537 | int iInt = chosen[i]; |
---|
538 | if (iInt>=0) { |
---|
539 | //std::cout<<"Strong branching on " |
---|
540 | // <<i<<""<<iInt<<" down "<<downMovement[i] |
---|
541 | // <<" up "<<upMovement[i] |
---|
542 | // <<" value "<<solutionValue[i] |
---|
543 | // <<std::endl; |
---|
544 | bool better = false; |
---|
545 | if (min(upMovement[i],downMovement[i])>best) { |
---|
546 | // smaller is better |
---|
547 | better=true; |
---|
548 | } else if (min(upMovement[i],downMovement[i])>best-1.0e-5) { |
---|
549 | if (max(upMovement[i],downMovement[i])>best2+1.0e-5) { |
---|
550 | // smaller is about same, but larger is better |
---|
551 | better=true; |
---|
552 | } |
---|
553 | } |
---|
554 | if (better) { |
---|
555 | best = min(upMovement[i],downMovement[i]); |
---|
556 | best2 = max(upMovement[i],downMovement[i]); |
---|
557 | variable_ = iInt; |
---|
558 | double value = solutionValue[i]; |
---|
559 | value = max(value,(double) lower_[variable_]); |
---|
560 | value = min(value,(double) upper_[variable_]); |
---|
561 | value_=value; |
---|
562 | if (upMovement[i]<=downMovement[i]) |
---|
563 | way_ = WAY_UP_DOWN__NOTHING_DONE; // up |
---|
564 | else |
---|
565 | way_ = WAY_DOWN_UP__NOTHING_DONE; // down |
---|
566 | } |
---|
567 | } |
---|
568 | } |
---|
569 | } |
---|
570 | // Delete the snapshot |
---|
571 | model.unmarkHotStart(); |
---|
572 | } |
---|
573 | ////// End of Strong branching |
---|
574 | #else |
---|
575 | variable_=-1; |
---|
576 | // This has hard coded integer tolerance |
---|
577 | double mostAway=INTEGER_TOLERANCE; |
---|
578 | int numberAway=0; |
---|
579 | for (i=0;i<numberIntegers;i++) { |
---|
580 | int iColumn = integer[i]; |
---|
581 | lower_[i]=(int)lower[iColumn]; |
---|
582 | upper_[i]=(int)upper[iColumn]; |
---|
583 | double value = solution[iColumn]; |
---|
584 | value = max(value,(double) lower_[i]); |
---|
585 | value = min(value,(double) upper_[i]); |
---|
586 | double nearest = floor(value+0.5); |
---|
587 | if (fabs(value-nearest)>INTEGER_TOLERANCE) |
---|
588 | numberAway++; |
---|
589 | if (fabs(value-nearest)>mostAway) { |
---|
590 | mostAway=fabs(value-nearest); |
---|
591 | variable_=i; |
---|
592 | value_=value; |
---|
593 | if (value<=nearest) |
---|
594 | way_ = WAY_UP_DOWN__NOTHING_DONE; // up |
---|
595 | else |
---|
596 | way_ = WAY_DOWN_UP__NOTHING_DONE; // down |
---|
597 | } |
---|
598 | } |
---|
599 | #endif |
---|
600 | } |
---|
601 | |
---|
602 | void |
---|
603 | DBNodeSimple::gutsOfCopy(const DBNodeSimple & rhs) |
---|
604 | { |
---|
605 | index_=rhs.index_; |
---|
606 | node_id_=rhs.node_id_; |
---|
607 | basis_ = rhs.basis_ ? rhs.basis_->clone() : NULL; |
---|
608 | objectiveValue_=rhs.objectiveValue_; |
---|
609 | variable_=rhs.variable_; |
---|
610 | way_=rhs.way_; |
---|
611 | numberIntegers_=rhs.numberIntegers_; |
---|
612 | value_=rhs.value_; |
---|
613 | parent_ = rhs.parent_; |
---|
614 | child_down_ = rhs.child_down_; |
---|
615 | child_up_ = rhs.child_up_; |
---|
616 | strong_branching_fixed_vars_ = rhs.strong_branching_fixed_vars_; |
---|
617 | reduced_cost_fixed_vars_ = rhs.reduced_cost_fixed_vars_; |
---|
618 | previous_ = rhs.previous_; |
---|
619 | next_ = rhs.next_; |
---|
620 | lower_=NULL; |
---|
621 | upper_=NULL; |
---|
622 | if (rhs.lower_!=NULL) { |
---|
623 | lower_ = new int [numberIntegers_]; |
---|
624 | upper_ = new int [numberIntegers_]; |
---|
625 | assert (upper_!=NULL); |
---|
626 | memcpy(lower_,rhs.lower_,numberIntegers_*sizeof(int)); |
---|
627 | memcpy(upper_,rhs.upper_,numberIntegers_*sizeof(int)); |
---|
628 | } |
---|
629 | } |
---|
630 | |
---|
631 | DBNodeSimple::DBNodeSimple(const DBNodeSimple & rhs) |
---|
632 | { |
---|
633 | gutsOfCopy(rhs); |
---|
634 | } |
---|
635 | |
---|
636 | DBNodeSimple & |
---|
637 | DBNodeSimple::operator=(const DBNodeSimple & rhs) |
---|
638 | { |
---|
639 | if (this != &rhs) { |
---|
640 | // LL: in original code basis/lower/upper was left alone if rhs did not |
---|
641 | // have them. Was that intentional? |
---|
642 | gutsOfDestructor(); |
---|
643 | gutsOfCopy(rhs); |
---|
644 | } |
---|
645 | return *this; |
---|
646 | } |
---|
647 | |
---|
648 | |
---|
649 | DBNodeSimple::~DBNodeSimple () |
---|
650 | { |
---|
651 | gutsOfDestructor(); |
---|
652 | } |
---|
653 | // Work of destructor |
---|
654 | void |
---|
655 | DBNodeSimple::gutsOfDestructor() |
---|
656 | { |
---|
657 | delete [] lower_; |
---|
658 | delete [] upper_; |
---|
659 | delete basis_; |
---|
660 | lower_ = NULL; |
---|
661 | upper_ = NULL; |
---|
662 | basis_ = NULL; |
---|
663 | objectiveValue_ = COIN_DBL_MAX; |
---|
664 | } |
---|
665 | // Extension - true if other extension of this |
---|
666 | bool |
---|
667 | DBNodeSimple::extension(const DBNodeSimple & other, |
---|
668 | const double * originalLower, |
---|
669 | const double * originalUpper) const |
---|
670 | { |
---|
671 | bool ok=true; |
---|
672 | for (int i=0;i<numberIntegers_;i++) { |
---|
673 | if (upper_[i]<originalUpper[i]|| |
---|
674 | lower_[i]>originalLower[i]) { |
---|
675 | if (other.upper_[i]>upper_[i]|| |
---|
676 | other.lower_[i]<lower_[i]) { |
---|
677 | ok=false; |
---|
678 | break; |
---|
679 | } |
---|
680 | } |
---|
681 | } |
---|
682 | return ok; |
---|
683 | } |
---|
684 | |
---|
685 | #include <vector> |
---|
686 | |
---|
687 | // Must code up by hand |
---|
688 | class DBVectorNode { |
---|
689 | |
---|
690 | public: |
---|
691 | |
---|
692 | // Default Constructor |
---|
693 | DBVectorNode (); |
---|
694 | |
---|
695 | // Copy constructor |
---|
696 | DBVectorNode ( const DBVectorNode &); |
---|
697 | |
---|
698 | // Assignment operator |
---|
699 | DBVectorNode & operator=( const DBVectorNode& rhs); |
---|
700 | |
---|
701 | // Destructor |
---|
702 | ~DBVectorNode (); |
---|
703 | // Size |
---|
704 | inline int size() const |
---|
705 | { return size_-sizeDeferred_;} |
---|
706 | // Push |
---|
707 | void push_back(DBNodeSimple & node); // the node_id_ of node will change |
---|
708 | /* Remove a single node from the tree and adjust the previos_/next_/first_ |
---|
709 | etc fields. Does NOT update child/parent relationship */ |
---|
710 | void removeNode(int n); |
---|
711 | /* Remove the subtree rooted at node n. properly adjusts previos_/next_ etc |
---|
712 | fields. Does NOT update child/parent relationships. */ |
---|
713 | void removeSubTree(int n); |
---|
714 | // Last one in (or other criterion) |
---|
715 | DBNodeSimple back() const; |
---|
716 | // Get rid of last one |
---|
717 | void pop_back(); |
---|
718 | // Works out best one |
---|
719 | int best() const; |
---|
720 | // Rearranges the tree |
---|
721 | void moveNodeUp(const int* which, |
---|
722 | OsiSolverInterface& model, DBNodeSimple & node); |
---|
723 | // Fix the bounds in the descendants of subroot |
---|
724 | void adjustBounds(int subroot, int brvar, int brlb, int brub); |
---|
725 | |
---|
726 | // Check that the bounds correspond to the branching decisions... |
---|
727 | void checkTree() const; |
---|
728 | void checkNode(int node) const; |
---|
729 | |
---|
730 | // Public data |
---|
731 | // Maximum size |
---|
732 | int maximumSize_; |
---|
733 | // Current size |
---|
734 | int size_; |
---|
735 | // Number still hanging around |
---|
736 | int sizeDeferred_; |
---|
737 | // First spare |
---|
738 | int firstSpare_; |
---|
739 | // First |
---|
740 | int first_; |
---|
741 | // Last |
---|
742 | int last_; |
---|
743 | // Chosen one |
---|
744 | mutable int chosen_; |
---|
745 | // Nodes |
---|
746 | DBNodeSimple * nodes_; |
---|
747 | }; |
---|
748 | |
---|
749 | |
---|
750 | DBVectorNode::DBVectorNode() : |
---|
751 | maximumSize_(10), |
---|
752 | size_(0), |
---|
753 | sizeDeferred_(0), |
---|
754 | firstSpare_(0), |
---|
755 | first_(-1), |
---|
756 | last_(-1) |
---|
757 | { |
---|
758 | nodes_ = new DBNodeSimple[maximumSize_]; |
---|
759 | for (int i=0;i<maximumSize_;i++) { |
---|
760 | nodes_[i].previous_=i-1; |
---|
761 | nodes_[i].next_=i+1; |
---|
762 | } |
---|
763 | } |
---|
764 | |
---|
765 | DBVectorNode::DBVectorNode(const DBVectorNode & rhs) |
---|
766 | { |
---|
767 | maximumSize_ = rhs.maximumSize_; |
---|
768 | size_ = rhs.size_; |
---|
769 | sizeDeferred_ = rhs.sizeDeferred_; |
---|
770 | firstSpare_ = rhs.firstSpare_; |
---|
771 | first_ = rhs.first_; |
---|
772 | last_ = rhs.last_; |
---|
773 | nodes_ = new DBNodeSimple[maximumSize_]; |
---|
774 | for (int i=0;i<maximumSize_;i++) { |
---|
775 | nodes_[i] = rhs.nodes_[i]; |
---|
776 | } |
---|
777 | } |
---|
778 | |
---|
779 | DBVectorNode & |
---|
780 | DBVectorNode::operator=(const DBVectorNode & rhs) |
---|
781 | { |
---|
782 | if (this != &rhs) { |
---|
783 | delete [] nodes_; |
---|
784 | maximumSize_ = rhs.maximumSize_; |
---|
785 | size_ = rhs.size_; |
---|
786 | sizeDeferred_ = rhs.sizeDeferred_; |
---|
787 | firstSpare_ = rhs.firstSpare_; |
---|
788 | first_ = rhs.first_; |
---|
789 | last_ = rhs.last_; |
---|
790 | nodes_ = new DBNodeSimple[maximumSize_]; |
---|
791 | for (int i=0;i<maximumSize_;i++) { |
---|
792 | nodes_[i] = rhs.nodes_[i]; |
---|
793 | } |
---|
794 | } |
---|
795 | return *this; |
---|
796 | } |
---|
797 | |
---|
798 | |
---|
799 | DBVectorNode::~DBVectorNode () |
---|
800 | { |
---|
801 | delete [] nodes_; |
---|
802 | } |
---|
803 | // Push |
---|
804 | void |
---|
805 | DBVectorNode::push_back(DBNodeSimple & node) |
---|
806 | { |
---|
807 | if (size_==maximumSize_) { |
---|
808 | assert (firstSpare_==size_); |
---|
809 | maximumSize_ = (maximumSize_*3)+10; |
---|
810 | DBNodeSimple * temp = new DBNodeSimple[maximumSize_]; |
---|
811 | int i; |
---|
812 | for (i=0;i<size_;i++) { |
---|
813 | temp[i]=nodes_[i]; |
---|
814 | } |
---|
815 | delete [] nodes_; |
---|
816 | nodes_ = temp; |
---|
817 | //firstSpare_=size_; |
---|
818 | int last = -1; |
---|
819 | for ( i=size_;i<maximumSize_;i++) { |
---|
820 | nodes_[i].previous_=last; |
---|
821 | nodes_[i].next_=i+1; |
---|
822 | last = i; |
---|
823 | } |
---|
824 | } |
---|
825 | assert (firstSpare_<maximumSize_); |
---|
826 | assert (nodes_[firstSpare_].previous_<0); |
---|
827 | int next = nodes_[firstSpare_].next_; |
---|
828 | nodes_[firstSpare_]=node; |
---|
829 | nodes_[firstSpare_].node_id_ = firstSpare_; |
---|
830 | if (last_>=0) { |
---|
831 | assert (nodes_[last_].next_==-1); |
---|
832 | nodes_[last_].next_=firstSpare_; |
---|
833 | } |
---|
834 | nodes_[firstSpare_].previous_=last_; |
---|
835 | nodes_[firstSpare_].next_=-1; |
---|
836 | if (last_==-1) { |
---|
837 | assert (first_==-1); |
---|
838 | first_ = firstSpare_; |
---|
839 | } |
---|
840 | last_=firstSpare_; |
---|
841 | if (next>=0&&next<maximumSize_) { |
---|
842 | firstSpare_ = next; |
---|
843 | nodes_[firstSpare_].previous_=-1; |
---|
844 | } else { |
---|
845 | firstSpare_=maximumSize_; |
---|
846 | } |
---|
847 | chosen_ = -1; |
---|
848 | //best(); |
---|
849 | size_++; |
---|
850 | if (node.bothChildDone()) |
---|
851 | sizeDeferred_++; |
---|
852 | } |
---|
853 | // Works out best one |
---|
854 | int |
---|
855 | DBVectorNode::best() const |
---|
856 | { |
---|
857 | // can modify |
---|
858 | chosen_=-1; |
---|
859 | if (chosen_<0) { |
---|
860 | chosen_=last_; |
---|
861 | while (nodes_[chosen_].bothChildDone()) { |
---|
862 | chosen_ = nodes_[chosen_].previous_; |
---|
863 | assert (chosen_>=0); |
---|
864 | } |
---|
865 | } |
---|
866 | return chosen_; |
---|
867 | } |
---|
868 | // Last one in (or other criterion) |
---|
869 | DBNodeSimple |
---|
870 | DBVectorNode::back() const |
---|
871 | { |
---|
872 | assert (last_>=0); |
---|
873 | return nodes_[best()]; |
---|
874 | } |
---|
875 | |
---|
876 | void |
---|
877 | DBVectorNode::removeNode(int n) |
---|
878 | { |
---|
879 | DBNodeSimple& node = nodes_[n]; |
---|
880 | if (node.bothChildDone()) |
---|
881 | sizeDeferred_--; |
---|
882 | int previous = node.previous_; |
---|
883 | int next = node.next_; |
---|
884 | node.~DBNodeSimple(); |
---|
885 | if (previous >= 0) { |
---|
886 | nodes_[previous].next_=next; |
---|
887 | } else { |
---|
888 | first_ = next; |
---|
889 | } |
---|
890 | if (next >= 0) { |
---|
891 | nodes_[next].previous_ = previous; |
---|
892 | } else { |
---|
893 | last_ = previous; |
---|
894 | } |
---|
895 | node.previous_ = -1; |
---|
896 | if (firstSpare_ >= 0) { |
---|
897 | node.next_ = firstSpare_; |
---|
898 | } else { |
---|
899 | node.next_ = -1; |
---|
900 | } |
---|
901 | firstSpare_ = n; |
---|
902 | assert (size_>0); |
---|
903 | size_--; |
---|
904 | } |
---|
905 | |
---|
906 | void |
---|
907 | DBVectorNode::removeSubTree(int n) |
---|
908 | { |
---|
909 | if (nodes_[n].child_down_ >= 0) { |
---|
910 | removeSubTree(nodes_[n].child_down_); |
---|
911 | } |
---|
912 | if (nodes_[n].child_up_ >= 0) { |
---|
913 | removeSubTree(nodes_[n].child_up_); |
---|
914 | } |
---|
915 | removeNode(n); |
---|
916 | } |
---|
917 | |
---|
918 | // Get rid of last one |
---|
919 | void |
---|
920 | DBVectorNode::pop_back() |
---|
921 | { |
---|
922 | // Temporary until more sophisticated |
---|
923 | //assert (last_==chosen_); |
---|
924 | removeNode(chosen_); |
---|
925 | chosen_ = -1; |
---|
926 | } |
---|
927 | |
---|
928 | bool |
---|
929 | DBNodeSimple::canSwitchParentWithGrandparent(const int* which, |
---|
930 | LPresult& lpres, |
---|
931 | const int * original_lower, |
---|
932 | const int * original_upper, |
---|
933 | DBVectorNode & branchingTree) |
---|
934 | { |
---|
935 | /* |
---|
936 | The current node ('this') is really the parent (P) and the solution in the |
---|
937 | lpres represents the child. The grandparent (GP) is this.parent_. Let's have |
---|
938 | variable names respecting the truth. |
---|
939 | */ |
---|
940 | #if !defined(FUNNY_BRANCHING) |
---|
941 | return false; |
---|
942 | #endif |
---|
943 | |
---|
944 | const int parent_id = node_id_; |
---|
945 | const DBNodeSimple& parent = branchingTree.nodes_[parent_id]; |
---|
946 | const int grandparent_id = parent.parent_; |
---|
947 | |
---|
948 | if (grandparent_id == -1) { |
---|
949 | // can't flip root higher up... |
---|
950 | return false; |
---|
951 | } |
---|
952 | const DBNodeSimple& grandparent = branchingTree.nodes_[grandparent_id]; |
---|
953 | |
---|
954 | // THINK: these are not going to happen (hopefully), but if they do, what |
---|
955 | // should we do??? |
---|
956 | if (lpres.isAbandoned) { |
---|
957 | printf("WARNING: lpres.isAbandoned true!\n"); |
---|
958 | return false; |
---|
959 | } |
---|
960 | if (lpres.isProvenDualInfeasible) { |
---|
961 | printf("WARNING: lpres.isProvenDualInfeasible true!\n"); |
---|
962 | return false; |
---|
963 | } |
---|
964 | if (lpres.isPrimalObjectiveLimitReached) { |
---|
965 | printf("WARNING: lpres.isPrimalObjectiveLimitReached true!\n"); |
---|
966 | return false; |
---|
967 | } |
---|
968 | |
---|
969 | if (parent.strong_branching_fixed_vars_ || parent.reduced_cost_fixed_vars_ || |
---|
970 | grandparent.strong_branching_fixed_vars_ || |
---|
971 | grandparent.reduced_cost_fixed_vars_) { |
---|
972 | return false; |
---|
973 | } |
---|
974 | |
---|
975 | const double direction = lpres.getObjSense; |
---|
976 | |
---|
977 | const int GP_brvar = grandparent.variable_; |
---|
978 | const int GP_brvar_fullid = which[GP_brvar]; |
---|
979 | const bool parent_is_down_child = parent_id == grandparent.child_down_; |
---|
980 | |
---|
981 | if (lpres.isIterationLimitReached) { |
---|
982 | // THINK: should we do anything? |
---|
983 | return false; |
---|
984 | } |
---|
985 | |
---|
986 | // At this point the only flags of interest are isProvenOptimal, |
---|
987 | // isDualObjectiveLimitReached and isProvenPrimalInfeasible. |
---|
988 | // Note that isDualObjectiveLimitReached can't be true alone (search for |
---|
989 | // mustResolve). |
---|
990 | |
---|
991 | if (lpres.isProvenOptimal) { |
---|
992 | // May or may not be over the dual obj limit. If the obj val of parent is |
---|
993 | // higher than that of the grandparent then we allow switching. |
---|
994 | // NOTE: grandparent's obj val can;t be 1e100 since then parent would not |
---|
995 | // exist... |
---|
996 | return false; |
---|
997 | if (lpres.getObjValue > grandparent.objectiveValue_ + 1e-8) { |
---|
998 | double djValue = lpres.getReducedCost[GP_brvar_fullid]*direction; |
---|
999 | if (djValue > 1.0e-6) { |
---|
1000 | // wants to go down |
---|
1001 | if (parent_is_down_child) { |
---|
1002 | return true; |
---|
1003 | } |
---|
1004 | if (lpres.getColLower[GP_brvar_fullid] > std::ceil(grandparent.value_)) { |
---|
1005 | return true; |
---|
1006 | } |
---|
1007 | } else if (djValue < -1.0e-6) { |
---|
1008 | // wants to go up |
---|
1009 | if (! parent_is_down_child) { |
---|
1010 | return true; |
---|
1011 | } |
---|
1012 | if (lpres.getColUpper[GP_brvar_fullid] < std::floor(grandparent.value_)) { |
---|
1013 | return true; |
---|
1014 | } |
---|
1015 | } |
---|
1016 | } |
---|
1017 | return false; |
---|
1018 | } |
---|
1019 | |
---|
1020 | // Problem is really infeasible and has a dual ray. |
---|
1021 | assert(lpres.isProvenPrimalInfeasible); |
---|
1022 | const int greatgrandparent_id = grandparent.parent_; |
---|
1023 | const int x = GP_brvar_fullid; // for easier reference... we'll use s_x |
---|
1024 | |
---|
1025 | /* |
---|
1026 | Now we are going to check that if we relax the GP's branching decision |
---|
1027 | then the child's LP relaxation is still infeasible. The test is done by |
---|
1028 | checking whether the column (or its negative) of the GP's branching |
---|
1029 | variable cuts off the dual ray proving the infeasibility. |
---|
1030 | */ |
---|
1031 | |
---|
1032 | const double* colDualRay = lpres.colDualRay; |
---|
1033 | const double ub_without_GP = greatgrandparent_id >= 0 ? |
---|
1034 | branchingTree.nodes_[greatgrandparent_id].upper_[GP_brvar] : |
---|
1035 | original_upper[GP_brvar]; |
---|
1036 | const double lb_without_GP = greatgrandparent_id >= 0 ? |
---|
1037 | branchingTree.nodes_[greatgrandparent_id].lower_[GP_brvar] : |
---|
1038 | original_lower[GP_brvar]; |
---|
1039 | |
---|
1040 | if (parent_is_down_child) { // upper bound on x is forced |
---|
1041 | if (colDualRay[x] < 1e-8) { // if yA_x] is non-positive then the var is |
---|
1042 | // at its lower bound (or has 0 red cost) |
---|
1043 | // so same dual ray still proves infeas |
---|
1044 | return true; |
---|
1045 | } |
---|
1046 | const double max_ub_increase = lpres.yb_plus_rl_minus_su / colDualRay[x]; |
---|
1047 | assert(max_ub_increase >= 1e-8); |
---|
1048 | const bool willSwitch = |
---|
1049 | max_ub_increase > ub_without_GP - lpres.getColUpper[x] + 1e-8; |
---|
1050 | if (willSwitch) { |
---|
1051 | // The switch will happen so we might as well update the "infeasibility" |
---|
1052 | lpres.yb_plus_rl_minus_su -= colDualRay[x] * (ub_without_GP - lpres.getColUpper[x]); |
---|
1053 | } |
---|
1054 | return willSwitch; |
---|
1055 | } else { // lower bound on x is forced |
---|
1056 | if (colDualRay[x] > -1e-8) { // if yA_x is non-negative then the var is |
---|
1057 | // at its upper bound (or has 0 red cost) |
---|
1058 | // so same dual ray still proves infeas |
---|
1059 | return true; |
---|
1060 | } |
---|
1061 | const double max_lb_decrease = - lpres.yb_plus_rl_minus_su / colDualRay[x]; |
---|
1062 | assert(max_lb_decrease >= 1e-8); |
---|
1063 | const bool willSwitch = |
---|
1064 | max_lb_decrease > lpres.getColLower[x] - lb_without_GP + 1e-8; |
---|
1065 | if (willSwitch) { |
---|
1066 | // The switch will happen so we might as well update the "infeasibility" |
---|
1067 | lpres.yb_plus_rl_minus_su += colDualRay[x] * (lpres.getColLower[x] - lb_without_GP); |
---|
1068 | } |
---|
1069 | return willSwitch; |
---|
1070 | } |
---|
1071 | // THINK: the above is definitely good for minimization problems. Is it good |
---|
1072 | // for max? |
---|
1073 | |
---|
1074 | return true; // to placate some compilers |
---|
1075 | } |
---|
1076 | |
---|
1077 | void |
---|
1078 | DBVectorNode::adjustBounds(int subroot, int brvar, int brlb, int brub) |
---|
1079 | { |
---|
1080 | assert(subroot != -1); |
---|
1081 | DBNodeSimple& node = nodes_[subroot]; |
---|
1082 | // Take the intersection of brvar's bounds in node and those passed in |
---|
1083 | brub = CoinMin(brub, node.upper_[brvar]); |
---|
1084 | brlb = CoinMax(brlb, node.lower_[brvar]); |
---|
1085 | if (brub < brlb) { |
---|
1086 | // This node became infeasible. Get rid of it and of its descendants |
---|
1087 | const int parent_id = node.parent_; |
---|
1088 | removeSubTree(subroot); |
---|
1089 | if (nodes_[parent_id].child_down_ == subroot) { |
---|
1090 | nodes_[parent_id].child_down_ = -1; |
---|
1091 | } else { |
---|
1092 | nodes_[parent_id].child_up_ = -1; |
---|
1093 | } |
---|
1094 | return; |
---|
1095 | } |
---|
1096 | if (node.variable_ == brvar) { |
---|
1097 | if (node.value_ < brlb) { |
---|
1098 | // child_down_ just became infeasible. Just cut out the current node |
---|
1099 | // together with its child_down_ from the tree and hang child_up_ on the |
---|
1100 | // parent of this node. |
---|
1101 | if (node.child_down_ >= 0) { |
---|
1102 | removeSubTree(node.child_down_); |
---|
1103 | } |
---|
1104 | const int parent_id = node.parent_; |
---|
1105 | const int child_remains = node.child_up_; |
---|
1106 | if (nodes_[parent_id].child_down_ == subroot) { |
---|
1107 | nodes_[parent_id].child_down_ = child_remains; |
---|
1108 | } else { |
---|
1109 | nodes_[parent_id].child_up_ = child_remains; |
---|
1110 | } |
---|
1111 | removeNode(subroot); |
---|
1112 | if (child_remains >= 0) { |
---|
1113 | nodes_[child_remains].parent_ = parent_id; |
---|
1114 | adjustBounds(child_remains, brvar, brlb, brub); |
---|
1115 | } |
---|
1116 | return; |
---|
1117 | } |
---|
1118 | if (node.value_ > brub) { |
---|
1119 | // child_up_ just became infeasible. Just cut out the current node |
---|
1120 | // together with its child_down_ from the tree and hang child_down_ on |
---|
1121 | // the parent of this node. |
---|
1122 | if (node.child_up_ >= 0) { |
---|
1123 | removeSubTree(node.child_up_); |
---|
1124 | } |
---|
1125 | const int parent_id = node.parent_; |
---|
1126 | const int child_remains = node.child_down_; |
---|
1127 | if (nodes_[parent_id].child_down_ == subroot) { |
---|
1128 | nodes_[parent_id].child_down_ = child_remains; |
---|
1129 | } else { |
---|
1130 | nodes_[parent_id].child_up_ = child_remains; |
---|
1131 | } |
---|
1132 | removeNode(subroot); |
---|
1133 | if (child_remains >= 0) { |
---|
1134 | nodes_[child_remains].parent_ = parent_id; |
---|
1135 | adjustBounds(child_remains, brvar, brlb, brub); |
---|
1136 | } |
---|
1137 | return; |
---|
1138 | } |
---|
1139 | // Now brlb < node.value_ < brub (value_ is fractional) |
---|
1140 | } |
---|
1141 | node.lower_[brvar] = brlb; |
---|
1142 | node.upper_[brvar] = brub; |
---|
1143 | if (node.child_down_ >= 0) { |
---|
1144 | adjustBounds(node.child_down_, brvar, brlb, brub); |
---|
1145 | } |
---|
1146 | if (node.child_up_ >= 0) { |
---|
1147 | adjustBounds(node.child_up_, brvar, brlb, brub); |
---|
1148 | } |
---|
1149 | } |
---|
1150 | |
---|
1151 | void |
---|
1152 | DBVectorNode::moveNodeUp(const int* which, |
---|
1153 | OsiSolverInterface& model, DBNodeSimple & node) |
---|
1154 | { |
---|
1155 | /* |
---|
1156 | The current node ('this') is really the parent (P). The grandparent (GP) |
---|
1157 | is this.parent_. Let's have variable names respecting the truth. |
---|
1158 | */ |
---|
1159 | const int parent_id = node.node_id_; |
---|
1160 | DBNodeSimple& parent = nodes_[parent_id]; |
---|
1161 | const int grandparent_id = parent.parent_; |
---|
1162 | assert(grandparent_id != -1); |
---|
1163 | DBNodeSimple& grandparent = nodes_[grandparent_id]; |
---|
1164 | const int greatgrandparent_id = grandparent.parent_; |
---|
1165 | |
---|
1166 | const bool parent_is_down_child = parent_id == grandparent.child_down_; |
---|
1167 | |
---|
1168 | #if defined(DEBUG_DYNAMIC_BRANCHING) |
---|
1169 | if (dyn_debug >= 100) { |
---|
1170 | printf("entered moveNodeUp\n"); |
---|
1171 | printf("parent_id %d grandparent_id %d greatgrandparent_id %d\n", |
---|
1172 | parent_id, grandparent_id, greatgrandparent_id); |
---|
1173 | printf("parent.way_ %d\n", parent.way_); |
---|
1174 | } |
---|
1175 | #endif |
---|
1176 | |
---|
1177 | |
---|
1178 | // First hang the nodes where they belong. |
---|
1179 | parent.parent_ = greatgrandparent_id; |
---|
1180 | grandparent.parent_ = parent_id; |
---|
1181 | const bool down_child_stays_with_parent = parent.workingOnDownChild(); |
---|
1182 | int& child_to_move = |
---|
1183 | down_child_stays_with_parent ? parent.child_up_ : parent.child_down_; |
---|
1184 | const bool child_to_move_is_processed = parent.bothChildDone(); |
---|
1185 | |
---|
1186 | #if defined(DEBUG_DYNAMIC_BRANCHING) |
---|
1187 | if (dyn_debug >= 1000) { |
---|
1188 | printf("parent_is_down_child %d down_child_stays_with_parent %d child_to_move %d\n", parent_is_down_child, down_child_stays_with_parent, child_to_move); |
---|
1189 | } |
---|
1190 | #endif |
---|
1191 | |
---|
1192 | if (parent_is_down_child) { |
---|
1193 | grandparent.child_down_ = child_to_move; |
---|
1194 | } else { |
---|
1195 | grandparent.child_up_ = child_to_move; |
---|
1196 | } |
---|
1197 | if (child_to_move >= 0) { |
---|
1198 | nodes_[child_to_move].parent_ = grandparent_id; |
---|
1199 | } |
---|
1200 | child_to_move = grandparent_id; |
---|
1201 | if (greatgrandparent_id >= 0) { |
---|
1202 | DBNodeSimple& greatgrandparent = nodes_[greatgrandparent_id]; |
---|
1203 | if (grandparent_id == greatgrandparent.child_down_) { |
---|
1204 | greatgrandparent.child_down_ = parent_id; |
---|
1205 | } else { |
---|
1206 | greatgrandparent.child_up_ = parent_id; |
---|
1207 | } |
---|
1208 | } |
---|
1209 | |
---|
1210 | #if defined(DEBUG_DYNAMIC_BRANCHING) |
---|
1211 | if (dyn_debug >= 1000) { |
---|
1212 | printf("after exchange\n"); |
---|
1213 | printf("parent.parent_ %d parent.child_down_ %d parent.child_up_ %d\n", |
---|
1214 | parent.parent_, parent.child_down_, parent.child_up_); |
---|
1215 | printf("grandparent.parent_ %d grandparent.child_down_ %d grandparent.child_up_ %d\n", |
---|
1216 | grandparent.parent_, grandparent.child_down_, grandparent.child_up_); |
---|
1217 | if (greatgrandparent_id >= 0) { |
---|
1218 | DBNodeSimple& greatgrandparent = nodes_[greatgrandparent_id]; |
---|
1219 | printf("greatgrandparent.parent_ %d greatgrandparent.child_down_ %d greatgrandparent.child_up_ %d\n", |
---|
1220 | greatgrandparent.parent_, greatgrandparent.child_down_, greatgrandparent.child_up_); |
---|
1221 | } |
---|
1222 | printf("exiting moveNodeUp\n"); |
---|
1223 | } |
---|
1224 | #endif |
---|
1225 | |
---|
1226 | |
---|
1227 | // Now make sure way_ is set properly |
---|
1228 | if (down_child_stays_with_parent) { |
---|
1229 | if (!parent.bothChildDone()) { |
---|
1230 | parent.way_ = DBNodeSimple::WAY_UP_DOWN__BOTH_DONE; |
---|
1231 | sizeDeferred_++; |
---|
1232 | } |
---|
1233 | if (grandparent.bothChildDone()) { |
---|
1234 | if (! child_to_move_is_processed) { |
---|
1235 | grandparent.way_ = parent_is_down_child ? |
---|
1236 | DBNodeSimple::WAY_UP_DOWN__UP_DONE : |
---|
1237 | DBNodeSimple::WAY_DOWN_UP__DOWN_DONE; |
---|
1238 | sizeDeferred_--; |
---|
1239 | } |
---|
1240 | } else { // only parent is processed from the two children of grandparent |
---|
1241 | if (! child_to_move_is_processed) { |
---|
1242 | // remove grandparent, none of its children is processed now, why |
---|
1243 | // force its branching decision? |
---|
1244 | removeNode(grandparent_id); |
---|
1245 | parent.child_up_ = -1; |
---|
1246 | parent.way_ = DBNodeSimple::WAY_DOWN_UP__DOWN_DONE; |
---|
1247 | sizeDeferred_--; |
---|
1248 | return; // No bound changes are needed on the GO side as GP is |
---|
1249 | // deleted... |
---|
1250 | } else { |
---|
1251 | grandparent.way_ = parent_is_down_child ? |
---|
1252 | DBNodeSimple::WAY_DOWN_UP__DOWN_DONE : |
---|
1253 | DBNodeSimple::WAY_UP_DOWN__UP_DONE; |
---|
1254 | } |
---|
1255 | } |
---|
1256 | } else { |
---|
1257 | if (!parent.bothChildDone()) { |
---|
1258 | parent.way_ = DBNodeSimple::WAY_DOWN_UP__BOTH_DONE; |
---|
1259 | sizeDeferred_++; |
---|
1260 | } |
---|
1261 | if (grandparent.bothChildDone()) { |
---|
1262 | if (! child_to_move_is_processed) { |
---|
1263 | grandparent.way_ = parent_is_down_child ? |
---|
1264 | DBNodeSimple::WAY_UP_DOWN__UP_DONE : |
---|
1265 | DBNodeSimple::WAY_DOWN_UP__DOWN_DONE; |
---|
1266 | sizeDeferred_--; |
---|
1267 | } |
---|
1268 | } else { // only parent is processed from the two children of grandparent |
---|
1269 | if (! child_to_move_is_processed) { |
---|
1270 | removeNode(grandparent_id); |
---|
1271 | parent.child_down_ = -1; |
---|
1272 | parent.way_ = DBNodeSimple::WAY_UP_DOWN__UP_DONE; |
---|
1273 | sizeDeferred_--; |
---|
1274 | return; |
---|
1275 | } else { |
---|
1276 | grandparent.way_ = parent_is_down_child ? |
---|
1277 | DBNodeSimple::WAY_DOWN_UP__DOWN_DONE : |
---|
1278 | DBNodeSimple::WAY_UP_DOWN__UP_DONE; |
---|
1279 | } |
---|
1280 | } |
---|
1281 | } |
---|
1282 | |
---|
1283 | // Now modify bounds |
---|
1284 | |
---|
1285 | // First, get rid of GP's bound change of its branching variable in the |
---|
1286 | // bound list of P. Note: If greatgrandparent_id is < 0 then GP is the root |
---|
1287 | // so its bounds are the original bounds. |
---|
1288 | const int GP_brvar = grandparent.variable_; |
---|
1289 | if (parent_is_down_child) { |
---|
1290 | const int old_upper = grandparent.upper_[GP_brvar]; |
---|
1291 | parent.upper_[GP_brvar] = old_upper; |
---|
1292 | if ((GP_brvar != parent.variable_) || |
---|
1293 | (GP_brvar == parent.variable_ && !down_child_stays_with_parent)) { |
---|
1294 | model.setColUpper(which[GP_brvar], old_upper); |
---|
1295 | } |
---|
1296 | } else { |
---|
1297 | const int old_lower = grandparent.lower_[GP_brvar]; |
---|
1298 | parent.lower_[GP_brvar] = old_lower; |
---|
1299 | if ((GP_brvar != parent.variable_) || |
---|
1300 | (GP_brvar == parent.variable_ && down_child_stays_with_parent)) { |
---|
1301 | model.setColLower(which[GP_brvar], old_lower); |
---|
1302 | } |
---|
1303 | } |
---|
1304 | |
---|
1305 | // Now add the branching var bound change of P to GP and all of its |
---|
1306 | // descendant |
---|
1307 | if (down_child_stays_with_parent) { |
---|
1308 | adjustBounds(grandparent_id, parent.variable_, |
---|
1309 | (int)ceil(parent.value_), parent.upper_[parent.variable_]); |
---|
1310 | } else { |
---|
1311 | adjustBounds(grandparent_id, parent.variable_, |
---|
1312 | parent.lower_[parent.variable_], (int)floor(parent.value_)); |
---|
1313 | } |
---|
1314 | } |
---|
1315 | |
---|
1316 | void |
---|
1317 | DBVectorNode::checkNode(int node) const |
---|
1318 | { |
---|
1319 | if (node == -1) { |
---|
1320 | return; |
---|
1321 | } |
---|
1322 | const DBNodeSimple* n = nodes_ + node; |
---|
1323 | const DBNodeSimple* p = nodes_ + n->parent_; |
---|
1324 | for (int i = n->numberIntegers_-1; i >= 0; --i) { |
---|
1325 | if (i == p->variable_) { |
---|
1326 | if (node == p->child_down_) { |
---|
1327 | assert(p->lower_[i] <= n->lower_[i]); |
---|
1328 | assert((int)(floor(p->value_)) == n->upper_[i]); |
---|
1329 | } else { |
---|
1330 | assert((int)(ceil(p->value_)) == n->lower_[i]); |
---|
1331 | assert(p->upper_[i] >= n->upper_[i]); |
---|
1332 | } |
---|
1333 | } else { |
---|
1334 | assert(p->lower_[i] <= n->lower_[i]); |
---|
1335 | assert(p->upper_[i] >= n->upper_[i]); |
---|
1336 | } |
---|
1337 | } |
---|
1338 | checkNode(n->child_down_); |
---|
1339 | checkNode(n->child_up_); |
---|
1340 | } |
---|
1341 | |
---|
1342 | void |
---|
1343 | DBVectorNode::checkTree() const |
---|
1344 | { |
---|
1345 | // find the root |
---|
1346 | int root = first_; |
---|
1347 | while (true) { |
---|
1348 | if (nodes_[root].parent_ == -1) { |
---|
1349 | break; |
---|
1350 | } else { |
---|
1351 | root = nodes_[root].parent_; |
---|
1352 | } |
---|
1353 | } |
---|
1354 | checkNode(nodes_[root].child_down_); |
---|
1355 | checkNode(nodes_[root].child_up_); |
---|
1356 | } |
---|
1357 | |
---|
1358 | std::string getTree(DBVectorNode& branchingTree) |
---|
1359 | { |
---|
1360 | std::string tree; |
---|
1361 | char line[1000]; |
---|
1362 | for(int k=0; k<branchingTree.size_; k++) { |
---|
1363 | DBNodeSimple& node = branchingTree.nodes_[k]; |
---|
1364 | sprintf(line, "%d %d %d %d %f %d %d 0x%x %d %d\n", |
---|
1365 | k, node.node_id_, node.parent_, node.variable_, |
---|
1366 | node.value_, node.lower_[node.variable_], |
---|
1367 | node.upper_[node.variable_], node.way_, |
---|
1368 | node.child_down_, node.child_up_); |
---|
1369 | tree += line; |
---|
1370 | } |
---|
1371 | return tree; |
---|
1372 | } |
---|
1373 | |
---|
1374 | void printTree(const std::string& tree, int levels) |
---|
1375 | { |
---|
1376 | size_t pos = tree.size(); |
---|
1377 | for (int i = levels-1; i >= 0; --i) { |
---|
1378 | pos = tree.rfind('\n', pos-1); |
---|
1379 | if (pos == std::string::npos) { |
---|
1380 | pos = 0; |
---|
1381 | break; |
---|
1382 | } |
---|
1383 | } |
---|
1384 | printf("%s", tree.c_str() + (pos+1)); |
---|
1385 | } |
---|
1386 | |
---|
1387 | void printChain(DBVectorNode& branchingTree, int k) |
---|
1388 | { |
---|
1389 | while (k != -1) { |
---|
1390 | DBNodeSimple& node = branchingTree.nodes_[k]; |
---|
1391 | printf(" %d %d %d %d %f %d %d 0x%x %d %d %c %c\n", |
---|
1392 | k, node.node_id_, node.parent_, node.variable_, |
---|
1393 | node.value_, node.lower_[node.variable_], |
---|
1394 | node.upper_[node.variable_], node.way_, |
---|
1395 | node.child_down_, node.child_up_, |
---|
1396 | node.strong_branching_fixed_vars_ ? 'T' : 'F', |
---|
1397 | node.reduced_cost_fixed_vars_ ? 'T' : 'F'); |
---|
1398 | k = node.parent_; |
---|
1399 | } |
---|
1400 | } |
---|
1401 | |
---|
1402 | // Invoke solver's built-in enumeration algorithm |
---|
1403 | void |
---|
1404 | branchAndBound(OsiSolverInterface & model) { |
---|
1405 | double time1 = CoinCpuTime(); |
---|
1406 | |
---|
1407 | OsiClpSolverInterface* clp = |
---|
1408 | dynamic_cast<OsiClpSolverInterface*>(&model); |
---|
1409 | if(clp) { |
---|
1410 | ClpSimplex* modelPtr = clp->getModelPtr(); |
---|
1411 | modelPtr->setLogLevel(0); |
---|
1412 | } |
---|
1413 | |
---|
1414 | // solve LP |
---|
1415 | model.initialSolve(); |
---|
1416 | |
---|
1417 | if (model.isProvenOptimal()&&!model.isDualObjectiveLimitReached()) { |
---|
1418 | // Continuous is feasible - find integers |
---|
1419 | int numberIntegers=0; |
---|
1420 | int numberColumns = model.getNumCols(); |
---|
1421 | int iColumn; |
---|
1422 | int i; |
---|
1423 | for (iColumn=0;iColumn<numberColumns;iColumn++) { |
---|
1424 | if( model.isInteger(iColumn)) |
---|
1425 | numberIntegers++; |
---|
1426 | } |
---|
1427 | if (!numberIntegers) { |
---|
1428 | std::cout<<"No integer variables" |
---|
1429 | <<std::endl; |
---|
1430 | return; |
---|
1431 | } |
---|
1432 | int * which = new int[numberIntegers]; // which variables are integer |
---|
1433 | // original bounds |
---|
1434 | int * originalLower = new int[numberIntegers]; |
---|
1435 | int * originalUpper = new int[numberIntegers]; |
---|
1436 | int * relaxedLower = new int[numberIntegers]; |
---|
1437 | int * relaxedUpper = new int[numberIntegers]; |
---|
1438 | { |
---|
1439 | const double * lower = model.getColLower(); |
---|
1440 | const double * upper = model.getColUpper(); |
---|
1441 | numberIntegers=0; |
---|
1442 | for (iColumn=0;iColumn<numberColumns;iColumn++) { |
---|
1443 | if( model.isInteger(iColumn)) { |
---|
1444 | originalLower[numberIntegers]=(int) lower[iColumn]; |
---|
1445 | originalUpper[numberIntegers]=(int) upper[iColumn]; |
---|
1446 | which[numberIntegers++]=iColumn; |
---|
1447 | } |
---|
1448 | } |
---|
1449 | } |
---|
1450 | double direction = model.getObjSense(); |
---|
1451 | // empty tree |
---|
1452 | DBVectorNode branchingTree; |
---|
1453 | |
---|
1454 | // Add continuous to it; |
---|
1455 | DBNodeSimple rootNode(0,model,numberIntegers,which,model.getWarmStart()); |
---|
1456 | // something extra may have been fixed by strong branching |
---|
1457 | // if so go round again |
---|
1458 | while (rootNode.variable_==numberIntegers) { |
---|
1459 | model.resolve(); |
---|
1460 | rootNode = DBNodeSimple(0,model,numberIntegers,which,model.getWarmStart()); |
---|
1461 | } |
---|
1462 | if (rootNode.objectiveValue_<1.0e100) { |
---|
1463 | // push on stack |
---|
1464 | branchingTree.push_back(rootNode); |
---|
1465 | } |
---|
1466 | |
---|
1467 | // For printing totals |
---|
1468 | int numberIterations=0; |
---|
1469 | int numberNodes =0; |
---|
1470 | DBNodeSimple bestNode; |
---|
1471 | ////// Start main while of branch and bound |
---|
1472 | // while until nothing on stack |
---|
1473 | while (branchingTree.size()) { |
---|
1474 | #if defined(DEBUG_DYNAMIC_BRANCHING) |
---|
1475 | if (dyn_debug >= 1000) { |
---|
1476 | printf("branchingTree.size = %d %d\n",branchingTree.size(),branchingTree.size_); |
---|
1477 | printf("i node_id parent child_down child_up\n"); |
---|
1478 | for(int k=0; k<branchingTree.size_; k++) { |
---|
1479 | DBNodeSimple& node = branchingTree.nodes_[k]; |
---|
1480 | printf("%d %d %d %d %d\n",k, node.node_id_, node.parent_, |
---|
1481 | node.child_down_, node.child_up_); |
---|
1482 | } |
---|
1483 | } |
---|
1484 | #endif |
---|
1485 | // last node |
---|
1486 | DBNodeSimple node = branchingTree.back(); |
---|
1487 | int kNode = branchingTree.chosen_; |
---|
1488 | branchingTree.pop_back(); |
---|
1489 | #if defined(DEBUG_DYNAMIC_BRANCHING) |
---|
1490 | if (dyn_debug >= 1000) { |
---|
1491 | printf("Deleted current parent %d %d\n",branchingTree.size(),branchingTree.size_); |
---|
1492 | } |
---|
1493 | #endif |
---|
1494 | assert (! node.bothChildDone()); |
---|
1495 | numberNodes++; |
---|
1496 | if (node.variable_ < 0) { |
---|
1497 | // put back on tree and pretend both children are done. We want the |
---|
1498 | // whole tree all the time. |
---|
1499 | node.way_ = DBNodeSimple::WAY_UP_DOWN__BOTH_DONE; |
---|
1500 | branchingTree.push_back(node); |
---|
1501 | // Integer solution - save |
---|
1502 | bestNode = node; |
---|
1503 | // set cutoff (hard coded tolerance) |
---|
1504 | const double limit = (bestNode.objectiveValue_-1.0e-5)*direction; |
---|
1505 | model.setDblParam(OsiDualObjectiveLimit, limit); |
---|
1506 | std::cout<<"Integer solution of " |
---|
1507 | <<bestNode.objectiveValue_ |
---|
1508 | <<" found after "<<numberIterations |
---|
1509 | <<" iterations and "<<numberNodes<<" nodes" |
---|
1510 | <<std::endl; |
---|
1511 | #if defined(DEBUG_DYNAMIC_BRANCHING) |
---|
1512 | if (dyn_debug >= 1) { |
---|
1513 | branchingTree.checkTree(); |
---|
1514 | } |
---|
1515 | #endif |
---|
1516 | |
---|
1517 | } else { |
---|
1518 | // branch - do bounds |
---|
1519 | for (i=0;i<numberIntegers;i++) { |
---|
1520 | iColumn=which[i]; |
---|
1521 | model.setColBounds( iColumn,node.lower_[i],node.upper_[i]); |
---|
1522 | } |
---|
1523 | // move basis |
---|
1524 | model.setWarmStart(node.basis_); |
---|
1525 | // do branching variable |
---|
1526 | const bool down_branch = node.advanceWay(); |
---|
1527 | if (down_branch) { |
---|
1528 | model.setColUpper(which[node.variable_],floor(node.value_)); |
---|
1529 | } else { |
---|
1530 | model.setColLower(which[node.variable_],ceil(node.value_)); |
---|
1531 | } |
---|
1532 | // put back on tree anyway regardless whether any processing is left |
---|
1533 | // to be done. We want the whole tree all the time. |
---|
1534 | branchingTree.push_back(node); |
---|
1535 | #if defined(DEBUG_DYNAMIC_BRANCHING) |
---|
1536 | if (dyn_debug >= 1000) { |
---|
1537 | printf("Added current parent %d %d\n",branchingTree.size(),branchingTree.size_); |
---|
1538 | } |
---|
1539 | #endif |
---|
1540 | |
---|
1541 | // solve |
---|
1542 | model.resolve(); |
---|
1543 | CoinWarmStart * ws = model.getWarmStart(); |
---|
1544 | const CoinWarmStartBasis* wsb = |
---|
1545 | dynamic_cast<const CoinWarmStartBasis*>(ws); |
---|
1546 | assert (wsb!=NULL); // make sure not volume |
---|
1547 | numberIterations += model.getIterationCount(); |
---|
1548 | // fix on reduced costs |
---|
1549 | int nFixed0=0,nFixed1=0; |
---|
1550 | double cutoff; |
---|
1551 | model.getDblParam(OsiDualObjectiveLimit,cutoff); |
---|
1552 | double gap=(cutoff-model.getObjValue())*direction+1.0e-4; |
---|
1553 | // double |
---|
1554 | // gap=(cutoff-modelPtr_->objectiveValue())*direction+1.0e-4; |
---|
1555 | bool did_reduced_cost_fixing_for_child = false; |
---|
1556 | if (gap<1.0e10&&model.isProvenOptimal()&&!model.isDualObjectiveLimitReached()) { |
---|
1557 | const double * dj = model.getReducedCost(); |
---|
1558 | const double * lower = model.getColLower(); |
---|
1559 | const double * upper = model.getColUpper(); |
---|
1560 | for (i=0;i<numberIntegers;i++) { |
---|
1561 | iColumn=which[i]; |
---|
1562 | if (upper[iColumn]>lower[iColumn]) { |
---|
1563 | double djValue = dj[iColumn]*direction; |
---|
1564 | if (wsb->getStructStatus(iColumn)==CoinWarmStartBasis::atLowerBound&& |
---|
1565 | djValue>gap) { |
---|
1566 | nFixed0++; |
---|
1567 | model.setColUpper(iColumn,lower[iColumn]); |
---|
1568 | } else if (wsb->getStructStatus(iColumn)==CoinWarmStartBasis::atUpperBound&& |
---|
1569 | -djValue>gap) { |
---|
1570 | nFixed1++; |
---|
1571 | model.setColLower(iColumn,upper[iColumn]); |
---|
1572 | } |
---|
1573 | } |
---|
1574 | } |
---|
1575 | if (nFixed0+nFixed1) { |
---|
1576 | // printf("%d fixed to lower, %d fixed to upper\n", |
---|
1577 | // nFixed0,nFixed1); |
---|
1578 | did_reduced_cost_fixing_for_child = true; |
---|
1579 | } |
---|
1580 | } |
---|
1581 | if (model.isAbandoned()) { |
---|
1582 | // THINK: What the heck should we do??? |
---|
1583 | abort(); |
---|
1584 | } |
---|
1585 | if (model.isIterationLimitReached()) { |
---|
1586 | // maximum iterations - exit |
---|
1587 | std::cout<<"Exiting on maximum iterations\n"; |
---|
1588 | break; |
---|
1589 | } |
---|
1590 | |
---|
1591 | const double parentGap = (cutoff-node.objectiveValue_)*direction + 1.0e-4; |
---|
1592 | assert (parentGap >= 0); |
---|
1593 | const bool smallGap = false; // parentGap / fabs(cutoff) < 0.05; |
---|
1594 | |
---|
1595 | // We are not going to do any switching unless the gap is small |
---|
1596 | if (smallGap) { |
---|
1597 | // THINK: If goes over the dual obj limit then clp sets dual obj limit |
---|
1598 | // reached *AND* poven primal infeas. So we need to run it to |
---|
1599 | // completition to know whether it's really infeas or just over the |
---|
1600 | // bound. Something faster would be better... |
---|
1601 | const bool mustResolve = |
---|
1602 | model.isDualObjectiveLimitReached() && !model.isProvenOptimal(); |
---|
1603 | double oldlimit = 0; |
---|
1604 | if (mustResolve) { |
---|
1605 | model.getDblParam(OsiDualObjectiveLimit, oldlimit); |
---|
1606 | model.setDblParam(OsiDualObjectiveLimit, 1e100); |
---|
1607 | model.resolve(); |
---|
1608 | } |
---|
1609 | LPresult lpres(model); |
---|
1610 | if (mustResolve) { |
---|
1611 | lpres.isDualObjectiveLimitReached = true; |
---|
1612 | model.setDblParam(OsiDualObjectiveLimit, oldlimit); |
---|
1613 | } |
---|
1614 | bool canSwitch = |
---|
1615 | node.canSwitchParentWithGrandparent(which, lpres, originalLower, |
---|
1616 | originalUpper, branchingTree); |
---|
1617 | int cnt = 0; |
---|
1618 | |
---|
1619 | #if defined(DEBUG_DYNAMIC_BRANCHING) |
---|
1620 | if (dyn_debug >= 1) { |
---|
1621 | branchingTree.checkTree(); |
---|
1622 | } |
---|
1623 | #endif |
---|
1624 | |
---|
1625 | #if defined(DEBUG_DYNAMIC_BRANCHING) |
---|
1626 | std::string tree0; |
---|
1627 | if (dyn_debug >= 10) { |
---|
1628 | if (canSwitch) { |
---|
1629 | tree0 = getTree(branchingTree); |
---|
1630 | } |
---|
1631 | } |
---|
1632 | #endif |
---|
1633 | |
---|
1634 | while (canSwitch) { |
---|
1635 | branchingTree.moveNodeUp(which, model, node); |
---|
1636 | canSwitch = node.canSwitchParentWithGrandparent(which, lpres, |
---|
1637 | originalLower, |
---|
1638 | originalUpper, |
---|
1639 | branchingTree); |
---|
1640 | #if defined(DEBUG_DYNAMIC_BRANCHING) |
---|
1641 | if (dyn_debug >= 1) { |
---|
1642 | branchingTree.checkTree(); |
---|
1643 | } |
---|
1644 | #endif |
---|
1645 | ++cnt; |
---|
1646 | } |
---|
1647 | if (cnt > 0) { |
---|
1648 | printf("Before switch: opt: %c inf: %c dual_bd: %c\n", |
---|
1649 | lpres.isProvenOptimal ? 'T' : 'F', |
---|
1650 | lpres.isProvenPrimalInfeasible ? 'T' : 'F', |
---|
1651 | lpres.isDualObjectiveLimitReached ? 'T' : 'F'); |
---|
1652 | model.resolve(); |
---|
1653 | // This is horribly looking... Get rid of it when properly |
---|
1654 | // debugged... |
---|
1655 | #if 1 |
---|
1656 | const bool mustResolve = |
---|
1657 | model.isDualObjectiveLimitReached() && !model.isProvenOptimal(); |
---|
1658 | double oldlimit = 0; |
---|
1659 | |
---|
1660 | if (mustResolve) { |
---|
1661 | // THINK: Something faster would be better... |
---|
1662 | model.getDblParam(OsiDualObjectiveLimit, oldlimit); |
---|
1663 | model.setDblParam(OsiDualObjectiveLimit, 1e100); |
---|
1664 | model.resolve(); |
---|
1665 | } |
---|
1666 | LPresult lpres1(model); |
---|
1667 | if (mustResolve) { |
---|
1668 | lpres.isDualObjectiveLimitReached = true; |
---|
1669 | model.setDblParam(OsiDualObjectiveLimit, oldlimit); |
---|
1670 | } |
---|
1671 | printf("After resolve: opt: %c inf: %c dual_bd: %c\n", |
---|
1672 | lpres1.isProvenOptimal ? 'T' : 'F', |
---|
1673 | lpres1.isProvenPrimalInfeasible ? 'T' : 'F', |
---|
1674 | lpres1.isDualObjectiveLimitReached ? 'T' : 'F'); |
---|
1675 | assert(lpres.isAbandoned == model.isAbandoned()); |
---|
1676 | assert(lpres.isDualObjectiveLimitReached == model.isDualObjectiveLimitReached()); |
---|
1677 | assert(lpres.isDualObjectiveLimitReached || |
---|
1678 | (lpres.isProvenOptimal == model.isProvenOptimal())); |
---|
1679 | assert(lpres.isDualObjectiveLimitReached || |
---|
1680 | (lpres.isProvenPrimalInfeasible == model.isProvenPrimalInfeasible())); |
---|
1681 | assert(!lpres.isProvenOptimal || ! model.isProvenOptimal() || |
---|
1682 | (lpres.isProvenOptimal && model.isProvenOptimal() && |
---|
1683 | lpres.getObjValue == model.getObjValue())); |
---|
1684 | #endif |
---|
1685 | printf("Finished moving node %d up by %i levels.\n", node.node_id_, cnt); |
---|
1686 | } |
---|
1687 | |
---|
1688 | #if defined(DEBUG_DYNAMIC_BRANCHING) |
---|
1689 | if (dyn_debug >= 10) { |
---|
1690 | if (cnt > 0) { |
---|
1691 | std::string tree1 = getTree(branchingTree); |
---|
1692 | printf("=====================================\n"); |
---|
1693 | printf("It can move node %i up. way_: 0x%x brvar: %i\n", |
---|
1694 | node.node_id_, node.way_, node.variable_); |
---|
1695 | printTree(tree0, cnt+10); |
---|
1696 | printf("=====================================\n"); |
---|
1697 | printf("Finished moving the node up by %i levels.\n", cnt); |
---|
1698 | printTree(tree1, cnt+10); |
---|
1699 | printf("=====================================\n"); |
---|
1700 | } |
---|
1701 | } |
---|
1702 | #endif |
---|
1703 | } |
---|
1704 | if ((numberNodes%10)==0) |
---|
1705 | printf("%d nodes, tree size %d\n", |
---|
1706 | numberNodes,branchingTree.size()); |
---|
1707 | if (CoinCpuTime()-time1>3600.0) { |
---|
1708 | printf("stopping after 3600 seconds\n"); |
---|
1709 | exit(77); |
---|
1710 | } |
---|
1711 | DBNodeSimple newNode(numberNodes, model,numberIntegers,which,ws); |
---|
1712 | // something extra may have been fixed by strong branching |
---|
1713 | // if so go round again |
---|
1714 | while (newNode.variable_==numberIntegers) { |
---|
1715 | model.resolve(); |
---|
1716 | newNode = DBNodeSimple(numberNodes, model,numberIntegers,which,model.getWarmStart()); |
---|
1717 | newNode.strong_branching_fixed_vars_ = true; |
---|
1718 | } |
---|
1719 | newNode.reduced_cost_fixed_vars_ = did_reduced_cost_fixing_for_child; |
---|
1720 | if (newNode.objectiveValue_<1.0e100) { |
---|
1721 | newNode.parent_ = kNode; |
---|
1722 | // push on stack |
---|
1723 | branchingTree.push_back(newNode); |
---|
1724 | #if defined(DEBUG_DYNAMIC_BRANCHING) |
---|
1725 | if (dyn_debug >= 1000) { |
---|
1726 | printf("Added current child %d %d\n",branchingTree.size(),branchingTree.size_); |
---|
1727 | } |
---|
1728 | #endif |
---|
1729 | if (branchingTree.nodes_[kNode].workingOnDownChild()) { |
---|
1730 | branchingTree.nodes_[kNode].child_down_ = branchingTree.last_; |
---|
1731 | } else { |
---|
1732 | branchingTree.nodes_[kNode].child_up_ = branchingTree.last_; |
---|
1733 | } |
---|
1734 | if (newNode.variable_>=0) { |
---|
1735 | assert (fabs(newNode.value_-floor(newNode.value_+0.5))>1.0e-6); |
---|
1736 | } |
---|
1737 | #if 0 |
---|
1738 | else { |
---|
1739 | // integer solution - save |
---|
1740 | bestNode = node; |
---|
1741 | // set cutoff (hard coded tolerance) |
---|
1742 | model.setDblParam(OsiDualObjectiveLimit,(bestNode.objectiveValue_-1.0e-5)*direction); |
---|
1743 | std::cout<<"Integer solution of " |
---|
1744 | <<bestNode.objectiveValue_ |
---|
1745 | <<" found after "<<numberIterations |
---|
1746 | <<" iterations and "<<numberNodes<<" nodes" |
---|
1747 | <<std::endl; |
---|
1748 | } |
---|
1749 | #endif |
---|
1750 | } |
---|
1751 | #if defined(DEBUG_DYNAMIC_BRANCHING) |
---|
1752 | if (dyn_debug >= 1) { |
---|
1753 | branchingTree.checkTree(); |
---|
1754 | } |
---|
1755 | #endif |
---|
1756 | } |
---|
1757 | } |
---|
1758 | ////// End main while of branch and bound |
---|
1759 | std::cout<<"Search took " |
---|
1760 | <<numberIterations |
---|
1761 | <<" iterations and "<<numberNodes<<" nodes" |
---|
1762 | <<std::endl; |
---|
1763 | if (bestNode.numberIntegers_) { |
---|
1764 | // we have a solution restore |
---|
1765 | // do bounds |
---|
1766 | for (i=0;i<numberIntegers;i++) { |
---|
1767 | iColumn=which[i]; |
---|
1768 | model.setColBounds( iColumn,bestNode.lower_[i],bestNode.upper_[i]); |
---|
1769 | } |
---|
1770 | // move basis |
---|
1771 | model.setWarmStart(bestNode.basis_); |
---|
1772 | // set cutoff so will be good (hard coded tolerance) |
---|
1773 | model.setDblParam(OsiDualObjectiveLimit,(bestNode.objectiveValue_+1.0e-5)*direction); |
---|
1774 | model.resolve(); |
---|
1775 | } else { |
---|
1776 | OsiClpSolverInterface* clp = |
---|
1777 | dynamic_cast<OsiClpSolverInterface*>(&model); |
---|
1778 | if (clp) { |
---|
1779 | ClpSimplex* modelPtr_ = clp->getModelPtr(); |
---|
1780 | modelPtr_->setProblemStatus(1); |
---|
1781 | } |
---|
1782 | } |
---|
1783 | delete [] which; |
---|
1784 | delete [] originalLower; |
---|
1785 | delete [] originalUpper; |
---|
1786 | delete [] relaxedLower; |
---|
1787 | delete [] relaxedUpper; |
---|
1788 | } else { |
---|
1789 | std::cout<<"The LP relaxation is infeasible" |
---|
1790 | <<std::endl; |
---|
1791 | OsiClpSolverInterface* clp = |
---|
1792 | dynamic_cast<OsiClpSolverInterface*>(&model); |
---|
1793 | if (clp) { |
---|
1794 | ClpSimplex* modelPtr_ = clp->getModelPtr(); |
---|
1795 | modelPtr_->setProblemStatus(1); |
---|
1796 | } |
---|
1797 | //throw CoinError("The LP relaxation is infeasible or too expensive", |
---|
1798 | //"branchAndBound", "OsiClpSolverInterface"); |
---|
1799 | } |
---|
1800 | } |
---|
1801 | |
---|
1802 | |
---|
1803 | |
---|
1804 | int main(int argc, char* argv[]) |
---|
1805 | { |
---|
1806 | OsiClpSolverInterface model; |
---|
1807 | model.readMps(argv[1]); |
---|
1808 | branchAndBound(model); |
---|
1809 | return 0; |
---|
1810 | } |
---|