1 | /* |
---|
2 | Copyright (C) 2007, Lou Hafer, International Business Machines Corporation |
---|
3 | and others. All Rights Reserved. |
---|
4 | |
---|
5 | This file is part of cbc-generic. |
---|
6 | */ |
---|
7 | |
---|
8 | #include <iostream> |
---|
9 | |
---|
10 | #include "CoinTime.hpp" |
---|
11 | |
---|
12 | #include "OsiSolverInterface.hpp" |
---|
13 | #include "OsiChooseVariable.hpp" |
---|
14 | |
---|
15 | #include "CglPreProcess.hpp" |
---|
16 | |
---|
17 | #include "CbcModel.hpp" |
---|
18 | #include "CbcCutGenerator.hpp" |
---|
19 | #include "CbcBranchActual.hpp" |
---|
20 | #include "CbcStrategy.hpp" |
---|
21 | |
---|
22 | #include "CbcGenCtlBlk.hpp" |
---|
23 | #include "CbcGenParam.hpp" |
---|
24 | #include "CbcGenCbcParam.hpp" |
---|
25 | |
---|
26 | #define CBC_TRACK_SOLVERS 1 |
---|
27 | // #define COIN_CBC_VERBOSITY 5 |
---|
28 | |
---|
29 | /* |
---|
30 | The support functions for the main branch-and-cut action routine. |
---|
31 | */ |
---|
32 | |
---|
33 | namespace { |
---|
34 | |
---|
35 | char svnid[] = "$Id: CbcGenBaB.cpp 1173 2009-06-04 09:44:10Z forrest $" ; |
---|
36 | |
---|
37 | /* |
---|
38 | A hack to fix variables based on reduced cost prior to branch-and-cut. Note |
---|
39 | that we're *not* looking at the integrality gap here. Given the reduced costs |
---|
40 | of the root relaxation, we're simply placing a bet that variables with really |
---|
41 | unfavourable reduced costs that are at their most favourable bound in the |
---|
42 | root relaxation will never move from that bound. |
---|
43 | |
---|
44 | For the standard OsiSolverInterface, this requires a bit of effort as the |
---|
45 | solution and bounds arrays are const and the functions to change them have |
---|
46 | incompatible interfaces. |
---|
47 | */ |
---|
48 | |
---|
49 | void reducedCostHack (OsiSolverInterface *osi, double threshold) |
---|
50 | |
---|
51 | { |
---|
52 | int numCols = osi->getNumCols() ; |
---|
53 | int i ; |
---|
54 | const double *lower = osi->getColLower() ; |
---|
55 | const double *upper = osi->getColUpper() ; |
---|
56 | const double *solution = osi->getColSolution() ; |
---|
57 | const double *dj = osi->getReducedCost() ; |
---|
58 | /* |
---|
59 | First task: scan the columns looking for variables that are at their |
---|
60 | favourable bound and have reduced cost that exceeds the threshold. Remember |
---|
61 | the column index and the value. |
---|
62 | */ |
---|
63 | double *chgBnds = new double [numCols] ; |
---|
64 | int *chgCols = new int [numCols] ; |
---|
65 | |
---|
66 | int numFixed = 0 ; |
---|
67 | for (i = 0 ; i < numCols ; i++) { |
---|
68 | if (osi->isInteger(i)) { |
---|
69 | double value = solution[i] ; |
---|
70 | if (value < lower[i] + 1.0e-5 && dj[i] > threshold) { |
---|
71 | chgCols[numFixed] = i ; |
---|
72 | chgBnds[numFixed] = lower[i] ; |
---|
73 | numFixed++ ; |
---|
74 | } else if (value > upper[i] - 1.0e-5 && dj[i] < -threshold) { |
---|
75 | chgCols[numFixed] = i ; |
---|
76 | chgBnds[numFixed] = upper[i] ; |
---|
77 | numFixed++ ; |
---|
78 | } |
---|
79 | } |
---|
80 | } |
---|
81 | /* |
---|
82 | Second task: For variables that we want to fix, we need to: |
---|
83 | * Prepare an array with the new lower and upper bounds for variables that |
---|
84 | will be fixed. setColSetBounds requires an array with column indices and |
---|
85 | an array with new values for both bounds. |
---|
86 | * Set the correct value in a copy of the current solution. setColSolution |
---|
87 | requires a complete solution. |
---|
88 | */ |
---|
89 | if (numFixed > 0) { |
---|
90 | double *newSoln = CoinCopyOfArray(solution, numCols) ; |
---|
91 | double *newBnds = new double [2*numFixed] ; |
---|
92 | double *bndPtr = &newBnds[0] ; |
---|
93 | for (i = 0 ; i < numFixed ; i++) { |
---|
94 | int j = chgCols[i] ; |
---|
95 | double val = chgBnds[i] ; |
---|
96 | *bndPtr++ = val ; |
---|
97 | *bndPtr++ = val ; |
---|
98 | newSoln[j] = val ; |
---|
99 | } |
---|
100 | osi->setColSetBounds(&chgCols[0], &chgCols[numFixed], &newBnds[0]) ; |
---|
101 | osi->setColSolution(&newSoln[0]) ; |
---|
102 | |
---|
103 | std::cout |
---|
104 | << "Reduced cost fixing prior to B&C: " << numFixed |
---|
105 | << " columns fixed." << std::endl ; |
---|
106 | |
---|
107 | delete[] newSoln ; |
---|
108 | delete[] newBnds ; |
---|
109 | } |
---|
110 | |
---|
111 | delete[] chgBnds ; |
---|
112 | delete[] chgCols ; |
---|
113 | |
---|
114 | return ; |
---|
115 | } |
---|
116 | |
---|
117 | /* |
---|
118 | Helper routine to solve a continuous relaxation and print something |
---|
119 | intelligent when the result is other than optimal. Returns true if the |
---|
120 | result is optimal, false otherwise. |
---|
121 | */ |
---|
122 | |
---|
123 | bool solveRelaxation (CbcModel *model) |
---|
124 | |
---|
125 | { |
---|
126 | OsiSolverInterface *osi = model->solver() ; |
---|
127 | |
---|
128 | model->initialSolve() ; |
---|
129 | |
---|
130 | if (!(osi->isProvenOptimal())) { |
---|
131 | bool reason = false ; |
---|
132 | if (osi->isProvenPrimalInfeasible()) { |
---|
133 | std::cout |
---|
134 | << "Continuous relaxation is primal infeasible." << std::endl ; |
---|
135 | reason = true ; |
---|
136 | } |
---|
137 | if (osi->isProvenDualInfeasible()) { |
---|
138 | std::cout |
---|
139 | << "Continuous relaxation is dual infeasible." << std::endl ; |
---|
140 | reason = true ; |
---|
141 | } |
---|
142 | if (osi->isIterationLimitReached()) { |
---|
143 | std::cout |
---|
144 | << "Continuous solver reached iteration limit." << std::endl ; |
---|
145 | reason = true ; |
---|
146 | } |
---|
147 | if (osi->isAbandoned()) { |
---|
148 | std::cout |
---|
149 | << "Continuous solver abandoned the problem." << std::endl ; |
---|
150 | reason = true ; |
---|
151 | } |
---|
152 | if (reason == false) { |
---|
153 | std::cout |
---|
154 | << "Continuous solver failed for unknown reason." << std::endl ; |
---|
155 | } |
---|
156 | return (false) ; |
---|
157 | } |
---|
158 | |
---|
159 | return (true) ; |
---|
160 | } |
---|
161 | |
---|
162 | |
---|
163 | /* |
---|
164 | Helper routine to establish a priority vector. |
---|
165 | */ |
---|
166 | |
---|
167 | void setupPriorities (CbcModel *model, CbcGenCtlBlk::BPControl how) |
---|
168 | |
---|
169 | { |
---|
170 | int numCols = model->getNumCols() ; |
---|
171 | int *sort = new int[numCols] ; |
---|
172 | double *dsort = new double[numCols] ; |
---|
173 | int *priority = new int[numCols] ; |
---|
174 | const double *objective = model->getObjCoefficients() ; |
---|
175 | int iColumn ; |
---|
176 | int n = 0 ; |
---|
177 | bool priorityOK = true ; |
---|
178 | |
---|
179 | for (iColumn = 0 ; iColumn < numCols ; iColumn++) { |
---|
180 | if (model->isInteger(iColumn)) { |
---|
181 | sort[n] = n ; |
---|
182 | if (how == CbcGenCtlBlk::BPCost) { |
---|
183 | dsort[n++] = -objective[iColumn] ; |
---|
184 | } else if (how == CbcGenCtlBlk::BPOrder) { |
---|
185 | dsort[n++] = iColumn ; |
---|
186 | } else { |
---|
187 | std::cerr |
---|
188 | << "setupPriorities: Unrecognised priority specification." |
---|
189 | << std::endl ; |
---|
190 | priorityOK = false ; |
---|
191 | } |
---|
192 | } |
---|
193 | } |
---|
194 | |
---|
195 | if (priorityOK) { |
---|
196 | CoinSort_2(dsort, dsort + n, sort) ; |
---|
197 | |
---|
198 | int level = 0 ; |
---|
199 | double last = -1.0e100 ; |
---|
200 | for (int i = 0 ; i < n ; i++) { |
---|
201 | int iPut = sort[i] ; |
---|
202 | if (dsort[i] != last) { |
---|
203 | level++ ; |
---|
204 | last = dsort[i] ; |
---|
205 | } |
---|
206 | priority[iPut] = level ; |
---|
207 | } |
---|
208 | |
---|
209 | model->passInPriorities(priority, false) ; |
---|
210 | } |
---|
211 | |
---|
212 | delete [] priority ; |
---|
213 | delete [] sort ; |
---|
214 | delete [] dsort ; |
---|
215 | |
---|
216 | return ; |
---|
217 | } |
---|
218 | |
---|
219 | |
---|
220 | /* |
---|
221 | Helper routine to install a batch of heuristics. Each call to getXXXHeuristic |
---|
222 | will return a pointer to the heuristic object in gen iff the heuristic is |
---|
223 | enabled. |
---|
224 | */ |
---|
225 | |
---|
226 | void installHeuristics (CbcGenCtlBlk *ctlBlk, CbcModel *model) |
---|
227 | |
---|
228 | { |
---|
229 | CbcGenCtlBlk::CGControl action ; |
---|
230 | CbcHeuristic *gen ; |
---|
231 | CbcTreeLocal *localTree ; |
---|
232 | /* |
---|
233 | FPump goes first because it only works before there's a solution. |
---|
234 | */ |
---|
235 | action = ctlBlk->getFPump(gen, model) ; |
---|
236 | if (action != CbcGenCtlBlk::CGOff) { |
---|
237 | model->addHeuristic(gen, "FPump") ; |
---|
238 | } |
---|
239 | action = ctlBlk->getRounding(gen, model) ; |
---|
240 | if (action != CbcGenCtlBlk::CGOff) { |
---|
241 | model->addHeuristic(gen, "Rounding") ; |
---|
242 | } |
---|
243 | action = ctlBlk->getCombine(gen, model) ; |
---|
244 | if (action != CbcGenCtlBlk::CGOff) { |
---|
245 | model->addHeuristic(gen, "Combine") ; |
---|
246 | } |
---|
247 | action = ctlBlk->getGreedyCover(gen, model) ; |
---|
248 | if (action != CbcGenCtlBlk::CGOff) { |
---|
249 | model->addHeuristic(gen, "GCov") ; |
---|
250 | } |
---|
251 | action = ctlBlk->getGreedyEquality(gen, model) ; |
---|
252 | if (action != CbcGenCtlBlk::CGOff) { |
---|
253 | model->addHeuristic(gen, "GEq") ; |
---|
254 | } |
---|
255 | /* |
---|
256 | This one's a bit different. We acquire the local tree and install it in the |
---|
257 | model. |
---|
258 | */ |
---|
259 | action = ctlBlk->getTreeLocal(localTree, model) ; |
---|
260 | if (action != CbcGenCtlBlk::CGOff) { |
---|
261 | model->passInTreeHandler(*localTree) ; |
---|
262 | } |
---|
263 | |
---|
264 | return ; |
---|
265 | } |
---|
266 | |
---|
267 | |
---|
268 | /* |
---|
269 | Helper routine to install cut generators. |
---|
270 | |
---|
271 | I need to install the new lift-and-project generator (LandP). Also need to |
---|
272 | figure out stored cuts. |
---|
273 | */ |
---|
274 | |
---|
275 | void installCutGenerators (CbcGenCtlBlk *ctlBlk, CbcModel *model) |
---|
276 | |
---|
277 | { |
---|
278 | int switches[20] ; |
---|
279 | int genCnt = 0 ; |
---|
280 | CbcGenCtlBlk::CGControl action ; |
---|
281 | CglCutGenerator *gen ; |
---|
282 | |
---|
283 | /* |
---|
284 | The magic numbers for the howOften parameter that determines how often the |
---|
285 | generator is invoked. -100 is disabled, -99 is root only, -98 will stay |
---|
286 | active only so long as it generates cuts that improve the objective. A value |
---|
287 | 1 <= k <= 90 means the generator will be called every k nodes. If k is |
---|
288 | negative, then it can be switched off if unproductive. If k is positive, |
---|
289 | it'll carry on regardless. |
---|
290 | */ |
---|
291 | int howOften[CbcGenCtlBlk::CGMarker] ; |
---|
292 | howOften[CbcGenCtlBlk::CGOff] = -100 ; |
---|
293 | howOften[CbcGenCtlBlk::CGOn] = -1 ; |
---|
294 | howOften[CbcGenCtlBlk::CGRoot] = -99 ; |
---|
295 | howOften[CbcGenCtlBlk::CGIfMove] = -98 ; |
---|
296 | howOften[CbcGenCtlBlk::CGForceOn] = 1 ; |
---|
297 | howOften[CbcGenCtlBlk::CGForceBut] = 1 ; |
---|
298 | |
---|
299 | /* |
---|
300 | A negative value for rowCuts means that the specified actions happen only at |
---|
301 | the root. |
---|
302 | */ |
---|
303 | action = ctlBlk->getProbing(gen) ; |
---|
304 | if (action != CbcGenCtlBlk::CGOff) { |
---|
305 | if (action == CbcGenCtlBlk::CGForceBut) { |
---|
306 | CglProbing *probingGen = dynamic_cast<CglProbing *>(gen) ; |
---|
307 | probingGen->setRowCuts(-3) ; |
---|
308 | } |
---|
309 | model->addCutGenerator(gen, howOften[action], "Probing") ; |
---|
310 | switches[genCnt++] = 0 ; |
---|
311 | } |
---|
312 | action = ctlBlk->getGomory(gen) ; |
---|
313 | if (action != CbcGenCtlBlk::CGOff) { |
---|
314 | model->addCutGenerator(gen, howOften[action], "Gomory") ; |
---|
315 | switches[genCnt++] = -1 ; |
---|
316 | } |
---|
317 | action = ctlBlk->getKnapsack(gen) ; |
---|
318 | if (action != CbcGenCtlBlk::CGOff) { |
---|
319 | model->addCutGenerator(gen, howOften[action], "Knapsack") ; |
---|
320 | switches[genCnt++] = 0 ; |
---|
321 | } |
---|
322 | action = ctlBlk->getRedSplit(gen) ; |
---|
323 | if (action != CbcGenCtlBlk::CGOff) { |
---|
324 | model->addCutGenerator(gen, howOften[action], "RedSplit") ; |
---|
325 | switches[genCnt++] = 1 ; |
---|
326 | } |
---|
327 | action = ctlBlk->getClique(gen) ; |
---|
328 | if (action != CbcGenCtlBlk::CGOff) { |
---|
329 | model->addCutGenerator(gen, howOften[action], "Clique") ; |
---|
330 | switches[genCnt++] = 0 ; |
---|
331 | } |
---|
332 | action = ctlBlk->getMir(gen) ; |
---|
333 | if (action != CbcGenCtlBlk::CGOff) { |
---|
334 | model->addCutGenerator(gen, howOften[action], "MIR2") ; |
---|
335 | switches[genCnt++] = -1 ; |
---|
336 | } |
---|
337 | action = ctlBlk->getFlow(gen) ; |
---|
338 | if (action != CbcGenCtlBlk::CGOff) { |
---|
339 | model->addCutGenerator(gen, howOften[action], "Flow") ; |
---|
340 | switches[genCnt++] = 1 ; |
---|
341 | } |
---|
342 | action = ctlBlk->getTwomir(gen) ; |
---|
343 | if (action != CbcGenCtlBlk::CGOff) { |
---|
344 | model->addCutGenerator(gen, howOften[action], "2-MIR") ; |
---|
345 | switches[genCnt++] = 1 ; |
---|
346 | } |
---|
347 | /* |
---|
348 | Set control parameters on cut generators. cutDepth says `use this generator |
---|
349 | when (depth in tree) mod cutDepth == 0'. setSwitchOffIfLessThan says `switch |
---|
350 | this generator off if the number of cuts at the root is less than the given |
---|
351 | value'. Sort of. I need to document the magic numbers for howOften , etc. |
---|
352 | */ |
---|
353 | genCnt = model->numberCutGenerators() ; |
---|
354 | int iGen ; |
---|
355 | for (iGen = 0 ; iGen < genCnt ; iGen++) { |
---|
356 | CbcCutGenerator *generator = model->cutGenerator(iGen) ; |
---|
357 | int howOften = generator->howOften() ; |
---|
358 | if (howOften == -98 || howOften == -99) { |
---|
359 | generator->setSwitchOffIfLessThan(switches[iGen]) ; |
---|
360 | } |
---|
361 | generator->setTiming(true) ; |
---|
362 | int cutDepth = ctlBlk->getCutDepth() ; |
---|
363 | if (cutDepth >= 0) { |
---|
364 | generator->setWhatDepth(cutDepth) ; |
---|
365 | } |
---|
366 | } |
---|
367 | /* |
---|
368 | Now some additional control parameters that affect cut generation activity. |
---|
369 | |
---|
370 | Minimum drop is the minimum objective degradation required to continue with |
---|
371 | cut passes. We want at least .05 unless the objective is tiny, in which |
---|
372 | case we'll drop down to a floor of .0001. |
---|
373 | */ |
---|
374 | { |
---|
375 | double objFrac = fabs(model->getMinimizationObjValue()) * .001 + .0001 ; |
---|
376 | double minDrop = CoinMin(.05, objFrac) ; |
---|
377 | model->setMinimumDrop(minDrop) ; |
---|
378 | } |
---|
379 | /* |
---|
380 | Set the maximum number of rounds of cut generation at the root and at nodes |
---|
381 | in the tree. If the value is positive, cut generation will terminate early |
---|
382 | if the objective degradation doesn't meet the minimum drop requirement. If |
---|
383 | the value is negatie, minimum drop is not considered. |
---|
384 | |
---|
385 | At the root, for small problems, push for 100 passes (really we're betting |
---|
386 | that we'll stop because no cuts were generated). For medium size problems, |
---|
387 | the same but say we can quit if we're not achieving the minimum drop. For |
---|
388 | big problems, cut the number of rounds to 20. The user may have expressed |
---|
389 | an opinion; if so, it's already set. |
---|
390 | |
---|
391 | Once we're in the tree, aim for one pass per activation. |
---|
392 | */ |
---|
393 | if (ctlBlk->setByUser_[CbcCbcParam::CUTPASS] == false) { |
---|
394 | int numCols = model->getNumCols() ; |
---|
395 | if (numCols < 500) |
---|
396 | model->setMaximumCutPassesAtRoot(-100) ; |
---|
397 | else if (numCols < 5000) |
---|
398 | model->setMaximumCutPassesAtRoot(100) ; |
---|
399 | else |
---|
400 | model->setMaximumCutPassesAtRoot(20) ; |
---|
401 | } |
---|
402 | |
---|
403 | model->setMaximumCutPasses(1) ; |
---|
404 | |
---|
405 | return ; |
---|
406 | } |
---|
407 | |
---|
408 | /* |
---|
409 | Install `objects' (integers, SOS sets, etc.) in the OSI. Cribbed from |
---|
410 | CoinSolve 061216 and subjected to moderate rewriting. A substantial amount |
---|
411 | of code that's only relevant for AMPL has been deleted. We're only supporting |
---|
412 | OsiObjects in cbc-generic. |
---|
413 | */ |
---|
414 | |
---|
415 | void setupObjects (OsiSolverInterface *osi, |
---|
416 | bool didIPP, CglPreProcess *ippObj) |
---|
417 | |
---|
418 | { |
---|
419 | int numInts = osi->getNumIntegers() ; |
---|
420 | int numObjs = osi->numberObjects() ; |
---|
421 | /* |
---|
422 | Does this OSI have defined objects already? If not, we'd best define the |
---|
423 | basic integer objects. |
---|
424 | */ |
---|
425 | if (numInts == 0 || numObjs == 0) { |
---|
426 | osi->findIntegers(false) ; |
---|
427 | numInts = osi->getNumIntegers() ; |
---|
428 | numObjs = osi->numberObjects() ; |
---|
429 | } |
---|
430 | /* |
---|
431 | If we did preprocessing and discovered SOS sets, create SOS objects and |
---|
432 | install them in the OSI. The priority of SOS objects is set so that larger |
---|
433 | sets have higher (lower numeric value) priority. The priority of the |
---|
434 | original objects is reset to be lower than the priority of any SOS object. |
---|
435 | Since the SOS objects are copied into the OSI, we need to delete our |
---|
436 | originals once they've been installed. |
---|
437 | |
---|
438 | It's not clear to me that this is the right thing to do, particularly if |
---|
439 | the OSI comes equipped with complex objects. -- lh, 061216 -- |
---|
440 | */ |
---|
441 | if (didIPP && ippObj->numberSOS()) { |
---|
442 | OsiObject **oldObjs = osi->objects() ; |
---|
443 | int numCols = osi->getNumCols() ; |
---|
444 | |
---|
445 | for (int iObj = 0 ; iObj < numObjs ; iObj++) { |
---|
446 | oldObjs[iObj]->setPriority(numCols + 1) ; |
---|
447 | } |
---|
448 | |
---|
449 | int numSOS = ippObj->numberSOS() ; |
---|
450 | OsiObject **sosObjs = new OsiObject *[numSOS] ; |
---|
451 | const int *starts = ippObj->startSOS() ; |
---|
452 | const int *which = ippObj->whichSOS() ; |
---|
453 | const int *type = ippObj->typeSOS() ; |
---|
454 | const double *weight = ippObj->weightSOS() ; |
---|
455 | int iSOS ; |
---|
456 | for (iSOS = 0 ; iSOS < numSOS ; iSOS++) { |
---|
457 | int iStart = starts[iSOS] ; |
---|
458 | int sosLen = starts[iSOS+1] - iStart ; |
---|
459 | sosObjs[iSOS] = |
---|
460 | new OsiSOS(osi, sosLen, which + iStart, weight + iStart, type[iSOS]) ; |
---|
461 | sosObjs[iSOS]->setPriority(numCols - sosLen) ; |
---|
462 | } |
---|
463 | osi->addObjects(numSOS, sosObjs) ; |
---|
464 | |
---|
465 | for (iSOS = 0 ; iSOS < numSOS ; iSOS++) |
---|
466 | delete sosObjs[iSOS] ; |
---|
467 | delete [] sosObjs ; |
---|
468 | } |
---|
469 | |
---|
470 | return ; |
---|
471 | } |
---|
472 | |
---|
473 | } // end local namespace |
---|
474 | |
---|
475 | |
---|
476 | namespace CbcGenParamUtils { |
---|
477 | |
---|
478 | /* |
---|
479 | Run branch-and-cut. |
---|
480 | */ |
---|
481 | |
---|
482 | int doBaCParam (CoinParam *param) |
---|
483 | |
---|
484 | { |
---|
485 | assert (param != 0) ; |
---|
486 | CbcGenParam *genParam = dynamic_cast<CbcGenParam *>(param) ; |
---|
487 | assert (genParam != 0) ; |
---|
488 | CbcGenCtlBlk *ctlBlk = genParam->obj() ; |
---|
489 | assert (ctlBlk != 0) ; |
---|
490 | CbcModel *model = ctlBlk->model_ ; |
---|
491 | assert (model != 0) ; |
---|
492 | /* |
---|
493 | Setup to return nonfatal/fatal error (1/-1) by default. |
---|
494 | */ |
---|
495 | int retval ; |
---|
496 | if (CoinParamUtils::isInteractive()) { |
---|
497 | retval = 1 ; |
---|
498 | } else { |
---|
499 | retval = -1 ; |
---|
500 | } |
---|
501 | ctlBlk->setBaBStatus(CbcGenCtlBlk::BACAbandon, CbcGenCtlBlk::BACmInvalid, |
---|
502 | CbcGenCtlBlk::BACwNotStarted, false, 0) ; |
---|
503 | /* |
---|
504 | We ain't gonna do squat without a good model. |
---|
505 | */ |
---|
506 | if (!ctlBlk->goodModel_) { |
---|
507 | std::cout << "** Current model not valid!" << std::endl ; |
---|
508 | return (retval) ; |
---|
509 | } |
---|
510 | /* |
---|
511 | Start the clock ticking. |
---|
512 | */ |
---|
513 | double time1 = CoinCpuTime() ; |
---|
514 | double time2 ; |
---|
515 | /* |
---|
516 | Create a clone of the model which we can modify with impunity. Extract |
---|
517 | the underlying solver for convenient access. |
---|
518 | */ |
---|
519 | CbcModel babModel(*model) ; |
---|
520 | OsiSolverInterface *babSolver = babModel.solver() ; |
---|
521 | assert (babSolver != 0) ; |
---|
522 | # if CBC_TRACK_SOLVERS > 0 |
---|
523 | std::cout |
---|
524 | << "doBaCParam: initial babSolver is " |
---|
525 | << std::hex << babSolver << std::dec |
---|
526 | << ", log level " << babSolver->messageHandler()->logLevel() |
---|
527 | << "." << std::endl ; |
---|
528 | # endif |
---|
529 | /* |
---|
530 | Solve the root relaxation. Bail unless it solves to optimality. |
---|
531 | */ |
---|
532 | if (!solveRelaxation(&babModel)) { |
---|
533 | ctlBlk->setBaBStatus(&babModel, CbcGenCtlBlk::BACwBareRoot) ; |
---|
534 | return (0) ; |
---|
535 | } |
---|
536 | # if COIN_CBC_VERBOSITY > 0 |
---|
537 | std::cout |
---|
538 | << "doBaCParam: initial relaxation z = " |
---|
539 | << babSolver->getObjValue() << "." << std::endl ; |
---|
540 | # endif |
---|
541 | /* |
---|
542 | Are we up for fixing variables based on reduced cost alone? |
---|
543 | */ |
---|
544 | if (ctlBlk->djFix_.action_ == true) { |
---|
545 | reducedCostHack(babSolver, ctlBlk->djFix_.threshold_) ; |
---|
546 | } |
---|
547 | /* |
---|
548 | Time to consider preprocessing. We'll do a bit of setup before getting to |
---|
549 | the meat of the issue. |
---|
550 | |
---|
551 | preIppSolver will hold a clone of the unpreprocessed constraint system. |
---|
552 | We'll need it when we postprocess. ippSolver holds the preprocessed |
---|
553 | constraint system. Again, we clone it and give the clone to babModel for |
---|
554 | B&C. Presumably we need an unmodified copy of the preprocessed system to |
---|
555 | do postprocessing, but the copy itself is hidden inside the preprocess |
---|
556 | object. |
---|
557 | */ |
---|
558 | OsiSolverInterface *preIppSolver = 0 ; |
---|
559 | CglPreProcess ippObj ; |
---|
560 | bool didIPP = false ; |
---|
561 | |
---|
562 | int numberChanged = 0 ; |
---|
563 | int numberOriginalColumns = babSolver->getNumCols() ; |
---|
564 | CbcGenCtlBlk::IPPControl ippAction = ctlBlk->getIPPAction() ; |
---|
565 | |
---|
566 | if (!(ippAction == CbcGenCtlBlk::IPPOff || |
---|
567 | ippAction == CbcGenCtlBlk::IPPStrategy)) { |
---|
568 | double timeLeft = babModel.getMaximumSeconds() ; |
---|
569 | preIppSolver = babSolver->clone() ; |
---|
570 | OsiSolverInterface *ippSolver ; |
---|
571 | # if CBC_TRACK_SOLVERS > 0 |
---|
572 | std::cout |
---|
573 | << "doBaCParam: clone made prior to IPP is " |
---|
574 | << std::hex << preIppSolver << std::dec |
---|
575 | << ", log level " << preIppSolver->messageHandler()->logLevel() |
---|
576 | << "." << std::endl ; |
---|
577 | # endif |
---|
578 | |
---|
579 | preIppSolver->setHintParam(OsiDoInBranchAndCut, true, OsiHintDo) ; |
---|
580 | ippObj.messageHandler()->setLogLevel(babModel.logLevel()) ; |
---|
581 | |
---|
582 | CglProbing probingGen ; |
---|
583 | probingGen.setUsingObjective(true) ; |
---|
584 | probingGen.setMaxPass(3) ; |
---|
585 | probingGen.setMaxProbeRoot(preIppSolver->getNumCols()) ; |
---|
586 | probingGen.setMaxElements(100) ; |
---|
587 | probingGen.setMaxLookRoot(50) ; |
---|
588 | probingGen.setRowCuts(3) ; |
---|
589 | ippObj.addCutGenerator(&probingGen) ; |
---|
590 | /* |
---|
591 | For preProcessNonDefault, the 2nd parameter controls the conversion of |
---|
592 | clique and SOS constraints. 0 does nothing, -1 converts <= to ==, and |
---|
593 | 2 and 3 form SOS sets under strict and not-so-strict conditions, |
---|
594 | respectively. |
---|
595 | */ |
---|
596 | int convert = 0 ; |
---|
597 | if (ippAction == CbcGenCtlBlk::IPPEqual) { |
---|
598 | convert = -1 ; |
---|
599 | } else if (ippAction == CbcGenCtlBlk::IPPEqualAll) { |
---|
600 | convert = -2 ; |
---|
601 | } else if (ippAction == CbcGenCtlBlk::IPPSOS) { |
---|
602 | convert = 2 ; |
---|
603 | } else if (ippAction == CbcGenCtlBlk::IPPTrySOS) { |
---|
604 | convert = 3 ; |
---|
605 | } |
---|
606 | |
---|
607 | ippSolver = ippObj.preProcessNonDefault(*preIppSolver, convert, 10) ; |
---|
608 | # if CBC_TRACK_SOLVERS > 0 |
---|
609 | std::cout |
---|
610 | << "doBaCParam: solver returned from IPP is " |
---|
611 | << std::hex << ippSolver << std::dec ; |
---|
612 | if (ippSolver) { |
---|
613 | std::cout |
---|
614 | << ", log level " << ippSolver->messageHandler()->logLevel() ; |
---|
615 | } |
---|
616 | std::cout << "." << std::endl ; |
---|
617 | # endif |
---|
618 | /* |
---|
619 | ippSolver == 0 is success of a sort --- integer preprocess has found the |
---|
620 | problem to be infeasible or unbounded. Need to think about how to indicate |
---|
621 | status. |
---|
622 | */ |
---|
623 | if (!ippSolver) { |
---|
624 | std::cout |
---|
625 | << "Integer preprocess says infeasible or unbounded" << std::endl ; |
---|
626 | delete preIppSolver ; |
---|
627 | ctlBlk->setBaBStatus(&babModel, CbcGenCtlBlk::BACwIPP) ; |
---|
628 | return (0) ; |
---|
629 | } |
---|
630 | # if COIN_CBC_VERBOSITY > 0 |
---|
631 | else { |
---|
632 | std::cout |
---|
633 | << "After integer preprocessing, model has " |
---|
634 | << ippSolver->getNumRows() |
---|
635 | << " rows, " << ippSolver->getNumCols() << " columns, and " |
---|
636 | << ippSolver->getNumElements() << " elements." << std::endl ; |
---|
637 | } |
---|
638 | # endif |
---|
639 | |
---|
640 | preIppSolver->setHintParam(OsiDoInBranchAndCut, false, OsiHintDo) ; |
---|
641 | ippSolver->setHintParam(OsiDoInBranchAndCut, false, OsiHintDo) ; |
---|
642 | |
---|
643 | if (ippAction == CbcGenCtlBlk::IPPSave) { |
---|
644 | ippSolver->writeMps("presolved", "mps", 1.0) ; |
---|
645 | std::cout |
---|
646 | << "Integer preprocessed model written to `presolved.mps' " |
---|
647 | << "as minimisation problem." << std::endl ; |
---|
648 | } |
---|
649 | |
---|
650 | OsiSolverInterface *osiTmp = ippSolver->clone() ; |
---|
651 | babModel.assignSolver(osiTmp) ; |
---|
652 | babSolver = babModel.solver() ; |
---|
653 | # if CBC_TRACK_SOLVERS > 0 |
---|
654 | std::cout |
---|
655 | << "doBaCParam: clone of IPP solver passed to babModel is " |
---|
656 | << std::hex << babSolver << std::dec |
---|
657 | << ", log level " << babSolver->messageHandler()->logLevel() |
---|
658 | << "." << std::endl ; |
---|
659 | # endif |
---|
660 | if (!solveRelaxation(&babModel)) { |
---|
661 | delete preIppSolver ; |
---|
662 | ctlBlk->setBaBStatus(&babModel, CbcGenCtlBlk::BACwIPPRelax) ; |
---|
663 | return (0) ; |
---|
664 | } |
---|
665 | # if COIN_CBC_VERBOSITY > 0 |
---|
666 | std::cout |
---|
667 | << "doBaCParam: presolved relaxation z = " |
---|
668 | << babSolver->getObjValue() << "." << std::endl ; |
---|
669 | # endif |
---|
670 | babModel.setMaximumSeconds(timeLeft - (CoinCpuTime() - time1)) ; |
---|
671 | didIPP = true ; |
---|
672 | } |
---|
673 | /* |
---|
674 | At this point, babModel and babSolver hold the constraint system we'll use |
---|
675 | for B&C (either the original system or the preprocessed system) and we have |
---|
676 | a solution to the lp relaxation. |
---|
677 | |
---|
678 | If we're using the COSTSTRATEGY option, set up priorities here and pass |
---|
679 | them to the babModel. |
---|
680 | */ |
---|
681 | if (ctlBlk->priorityAction_ != CbcGenCtlBlk::BPOff) { |
---|
682 | setupPriorities(&babModel, ctlBlk->priorityAction_) ; |
---|
683 | } |
---|
684 | /* |
---|
685 | Install heuristics and cutting planes. |
---|
686 | */ |
---|
687 | installHeuristics(ctlBlk, &babModel) ; |
---|
688 | installCutGenerators(ctlBlk, &babModel) ; |
---|
689 | /* |
---|
690 | Set up status print frequency for babModel. |
---|
691 | */ |
---|
692 | if (babModel.getNumCols() > 2000 || babModel.getNumRows() > 1500 || |
---|
693 | babModel.messageHandler()->logLevel() > 1) |
---|
694 | babModel.setPrintFrequency(100) ; |
---|
695 | /* |
---|
696 | If we've read in a known good solution for debugging, activate the row cut |
---|
697 | debugger. |
---|
698 | */ |
---|
699 | if (ctlBlk->debugSol_.values_) { |
---|
700 | if (ctlBlk->debugSol_.numCols_ == babModel.getNumCols()) { |
---|
701 | babSolver->activateRowCutDebugger(ctlBlk->debugSol_.values_) ; |
---|
702 | } else { |
---|
703 | std::cout |
---|
704 | << "doBaCParam: debug file has incorrect number of columns." |
---|
705 | << std::endl ; |
---|
706 | } |
---|
707 | } |
---|
708 | /* |
---|
709 | Set ratio-based integrality gap, if specified by user. |
---|
710 | */ |
---|
711 | if (ctlBlk->setByUser_[CbcCbcParam::GAPRATIO] == true) { |
---|
712 | double obj = babSolver->getObjValue() ; |
---|
713 | double gapRatio = babModel.getDblParam(CbcModel::CbcAllowableFractionGap) ; |
---|
714 | double gap = gapRatio * (1.0e-5 + fabs(obj)) ; |
---|
715 | babModel.setAllowableGap(gap) ; |
---|
716 | std::cout |
---|
717 | << "doBaCParam: Continuous objective = " << obj |
---|
718 | << ", so allowable gap set to " << gap << std::endl ; |
---|
719 | } |
---|
720 | /* |
---|
721 | A bit of mystery code. As best I can figure, setSpecialOptions(2) suppresses |
---|
722 | the removal of warm start information when checkSolution runs an lp to check |
---|
723 | a solution. John's comment, ``probably faster to use a basis to get integer |
---|
724 | solutions'' makes some sense in this context. Didn't try to track down |
---|
725 | moreMipOptions just yet. |
---|
726 | */ |
---|
727 | babModel.setSpecialOptions(babModel.specialOptions() | 2) ; |
---|
728 | /* |
---|
729 | { int ndx = whichParam(MOREMIPOPTIONS,numberParameters,parameters) ; |
---|
730 | int moreMipOptions = parameters[ndx].intValue() ; |
---|
731 | if (moreMipOptions >= 0) |
---|
732 | { printf("more mip options %d\n",moreMipOptions); |
---|
733 | babModel.setSearchStrategy(moreMipOptions); } } |
---|
734 | */ |
---|
735 | /* |
---|
736 | Begin the final run-up to branch-and-cut. |
---|
737 | |
---|
738 | Make sure that objects are set up in the solver. It's possible that whoever |
---|
739 | loaded the model into the solver also set up objects. But it's also |
---|
740 | entirely likely that none exist to this point (and interesting to note that |
---|
741 | IPP doesn't need to know anything about objects). |
---|
742 | */ |
---|
743 | setupObjects(babSolver, didIPP, &ippObj) ; |
---|
744 | /* |
---|
745 | Set the branching method. We can't do this until we establish objects, |
---|
746 | because the constructor will set up arrays based on the number of objects, |
---|
747 | and there's no provision to set this information after creation. Arguably not |
---|
748 | good --- it'd be nice to set this in the prototype model that's cloned for |
---|
749 | this routine. In CoinSolve, shadowPriceMode is handled with the TESTOSI |
---|
750 | option. |
---|
751 | */ |
---|
752 | OsiChooseStrong strong(babSolver) ; |
---|
753 | strong.setNumberBeforeTrusted(babModel.numberBeforeTrust()) ; |
---|
754 | strong.setNumberStrong(babModel.numberStrong()) ; |
---|
755 | strong.setShadowPriceMode(ctlBlk->chooseStrong_.shadowPriceMode_) ; |
---|
756 | CbcBranchDefaultDecision decision ; |
---|
757 | decision.setChooseMethod(strong) ; |
---|
758 | babModel.setBranchingMethod(decision) ; |
---|
759 | /* |
---|
760 | Here I've deleted a huge block of code that deals with external priorities, |
---|
761 | branch direction, pseudocosts, and solution. (PRIORITYIN) Also a block of |
---|
762 | code that generates C++ code. |
---|
763 | */ |
---|
764 | /* |
---|
765 | Set up strategy for branch-and-cut. Note that the integer code supplied to |
---|
766 | setupPreProcessing is *not* compatible with the IPPAction enum. But at least |
---|
767 | it's documented. See desiredPreProcess_ in CbcStrategyDefault. `1' is |
---|
768 | accidentally equivalent to IPPOn. |
---|
769 | */ |
---|
770 | |
---|
771 | if (ippAction == CbcGenCtlBlk::IPPStrategy) { |
---|
772 | CbcStrategyDefault strategy(true, 5, 5) ; |
---|
773 | strategy.setupPreProcessing(1) ; |
---|
774 | babModel.setStrategy(strategy) ; |
---|
775 | } |
---|
776 | /* |
---|
777 | Yes! At long last, we're ready for the big call. Do branch and cut. In |
---|
778 | general, the solver used to return the solution will not be the solver we |
---|
779 | passed in, so reset babSolver here. |
---|
780 | */ |
---|
781 | int statistics = (ctlBlk->printOpt_ > 0) ? ctlBlk->printOpt_ : 0 ; |
---|
782 | # if CBC_TRACK_SOLVERS > 0 |
---|
783 | std::cout |
---|
784 | << "doBaCParam: solver at call to branchAndBound is " |
---|
785 | << std::hex << babModel.solver() << std::dec |
---|
786 | << ", log level " << babModel.solver()->messageHandler()->logLevel() |
---|
787 | << "." << std::endl ; |
---|
788 | # endif |
---|
789 | babModel.branchAndBound(statistics) ; |
---|
790 | babSolver = babModel.solver() ; |
---|
791 | # if CBC_TRACK_SOLVERS > 0 |
---|
792 | std::cout |
---|
793 | << "doBaCParam: solver at return from branchAndBound is " |
---|
794 | << std::hex << babModel.solver() << std::dec |
---|
795 | << ", log level " << babModel.solver()->messageHandler()->logLevel() |
---|
796 | << "." << std::endl ; |
---|
797 | # endif |
---|
798 | /* |
---|
799 | Write out solution to preprocessed model. |
---|
800 | */ |
---|
801 | if (ctlBlk->debugCreate_ == "createAfterPre" && |
---|
802 | babModel.bestSolution()) { |
---|
803 | CbcGenParamUtils::saveSolution(babSolver, "debug.file") ; |
---|
804 | } |
---|
805 | /* |
---|
806 | Print some information about branch-and-cut. |
---|
807 | */ |
---|
808 | # if COIN_CBC_VERBOSITY > 0 |
---|
809 | std::cout |
---|
810 | << "Cuts at root node changed objective from " |
---|
811 | << babModel.getContinuousObjective() |
---|
812 | << " to " << babModel.rootObjectiveAfterCuts() << std::endl ; |
---|
813 | |
---|
814 | for (int iGen = 0 ; iGen < babModel.numberCutGenerators() ; iGen++) { |
---|
815 | CbcCutGenerator *generator = babModel.cutGenerator(iGen) ; |
---|
816 | std::cout |
---|
817 | << generator->cutGeneratorName() << " was tried " |
---|
818 | << generator->numberTimesEntered() << " times and created " |
---|
819 | << generator->numberCutsInTotal() << " cuts of which " |
---|
820 | << generator->numberCutsActive() |
---|
821 | << " were active after adding rounds of cuts" ; |
---|
822 | if (generator->timing()) { |
---|
823 | std::cout << " ( " << generator->timeInCutGenerator() << " seconds)" ; |
---|
824 | } |
---|
825 | std::cout << std::endl ; |
---|
826 | } |
---|
827 | # endif |
---|
828 | |
---|
829 | time2 = CoinCpuTime(); |
---|
830 | ctlBlk->totalTime_ += time2 - time1; |
---|
831 | /* |
---|
832 | If we performed integer preprocessing, time to back it out. |
---|
833 | */ |
---|
834 | if (ippAction != CbcGenCtlBlk::IPPOff) { |
---|
835 | # if CBC_TRACK_SOLVERS > 0 |
---|
836 | std::cout |
---|
837 | << "doBaCParam: solver passed to IPP postprocess is " |
---|
838 | << std::hex << babSolver << std::dec << "." << std::endl ; |
---|
839 | # endif |
---|
840 | ippObj.postProcess(*babSolver); |
---|
841 | babModel.assignSolver(preIppSolver) ; |
---|
842 | babSolver = babModel.solver() ; |
---|
843 | # if CBC_TRACK_SOLVERS > 0 |
---|
844 | std::cout |
---|
845 | << "doBaCParam: solver in babModel after IPP postprocess is " |
---|
846 | << std::hex << babSolver << std::dec << "." << std::endl ; |
---|
847 | # endif |
---|
848 | } |
---|
849 | /* |
---|
850 | Write out postprocessed solution to debug file, if requested. |
---|
851 | */ |
---|
852 | if (ctlBlk->debugCreate_ == "create" && babModel.bestSolution()) { |
---|
853 | CbcGenParamUtils::saveSolution(babSolver, "debug.file") ; |
---|
854 | } |
---|
855 | /* |
---|
856 | If we have a good solution, detach the solver with the answer. Fill in the |
---|
857 | rest of the status information for the benefit of the wider world. |
---|
858 | */ |
---|
859 | bool keepAnswerSolver = false ; |
---|
860 | OsiSolverInterface *answerSolver = 0 ; |
---|
861 | if (babModel.bestSolution()) { |
---|
862 | babModel.setModelOwnsSolver(false) ; |
---|
863 | keepAnswerSolver = true ; |
---|
864 | answerSolver = babSolver ; |
---|
865 | } |
---|
866 | ctlBlk->setBaBStatus(&babModel, CbcGenCtlBlk::BACwBAC, |
---|
867 | keepAnswerSolver, answerSolver) ; |
---|
868 | /* |
---|
869 | And one last bit of information & statistics. |
---|
870 | */ |
---|
871 | ctlBlk->printBaBStatus() ; |
---|
872 | std::cout << " " ; |
---|
873 | if (keepAnswerSolver) { |
---|
874 | std::cout |
---|
875 | << "objective " << babModel.getObjValue() << "; " ; |
---|
876 | } |
---|
877 | std::cout |
---|
878 | << babModel.getNodeCount() << " nodes and " |
---|
879 | << babModel.getIterationCount() << " iterations - took " |
---|
880 | << time2 - time1 << " seconds" << std::endl ; |
---|
881 | |
---|
882 | return (0) ; |
---|
883 | } |
---|
884 | |
---|
885 | } // end namespace CbcGenParamutils |
---|
886 | |
---|