1 | /* $Id: CbcCutGenerator.cpp 2105 2015-01-05 13:11:11Z forrest $ */ |
---|
2 | // Copyright (C) 2003, International Business Machines |
---|
3 | // Corporation and others. All Rights Reserved. |
---|
4 | // This code is licensed under the terms of the Eclipse Public License (EPL). |
---|
5 | |
---|
6 | #if defined(_MSC_VER) |
---|
7 | // Turn off compiler warning about long names |
---|
8 | # pragma warning(disable:4786) |
---|
9 | #endif |
---|
10 | #include "CbcConfig.h" |
---|
11 | #include <cassert> |
---|
12 | #include <cstdlib> |
---|
13 | #include <cmath> |
---|
14 | #include <cfloat> |
---|
15 | |
---|
16 | #ifdef COIN_HAS_CLP |
---|
17 | #include "OsiClpSolverInterface.hpp" |
---|
18 | #else |
---|
19 | #include "OsiSolverInterface.hpp" |
---|
20 | #endif |
---|
21 | //#define CGL_DEBUG 1 |
---|
22 | #ifdef CGL_DEBUG |
---|
23 | #include "OsiRowCutDebugger.hpp" |
---|
24 | #endif |
---|
25 | #include "CbcModel.hpp" |
---|
26 | #include "CbcMessage.hpp" |
---|
27 | #include "CbcCutGenerator.hpp" |
---|
28 | #include "CbcBranchDynamic.hpp" |
---|
29 | #include "CglProbing.hpp" |
---|
30 | #include "CoinTime.hpp" |
---|
31 | |
---|
32 | // Default Constructor |
---|
33 | CbcCutGenerator::CbcCutGenerator () |
---|
34 | : timeInCutGenerator_(0.0), |
---|
35 | model_(NULL), |
---|
36 | generator_(NULL), |
---|
37 | generatorName_(NULL), |
---|
38 | whenCutGenerator_(-1), |
---|
39 | whenCutGeneratorInSub_(-100), |
---|
40 | switchOffIfLessThan_(0), |
---|
41 | depthCutGenerator_(-1), |
---|
42 | depthCutGeneratorInSub_(-1), |
---|
43 | inaccuracy_(0), |
---|
44 | numberTimes_(0), |
---|
45 | numberCuts_(0), |
---|
46 | numberElements_(0), |
---|
47 | numberColumnCuts_(0), |
---|
48 | numberCutsActive_(0), |
---|
49 | numberCutsAtRoot_(0), |
---|
50 | numberActiveCutsAtRoot_(0), |
---|
51 | numberShortCutsAtRoot_(0), |
---|
52 | switches_(1), |
---|
53 | maximumTries_(-1) |
---|
54 | { |
---|
55 | } |
---|
56 | // Normal constructor |
---|
57 | CbcCutGenerator::CbcCutGenerator(CbcModel * model, CglCutGenerator * generator, |
---|
58 | int howOften, const char * name, |
---|
59 | bool normal, bool atSolution, |
---|
60 | bool infeasible, int howOftenInSub, |
---|
61 | int whatDepth, int whatDepthInSub, |
---|
62 | int switchOffIfLessThan) |
---|
63 | : |
---|
64 | timeInCutGenerator_(0.0), |
---|
65 | depthCutGenerator_(whatDepth), |
---|
66 | depthCutGeneratorInSub_(whatDepthInSub), |
---|
67 | inaccuracy_(0), |
---|
68 | numberTimes_(0), |
---|
69 | numberCuts_(0), |
---|
70 | numberElements_(0), |
---|
71 | numberColumnCuts_(0), |
---|
72 | numberCutsActive_(0), |
---|
73 | numberCutsAtRoot_(0), |
---|
74 | numberActiveCutsAtRoot_(0), |
---|
75 | numberShortCutsAtRoot_(0), |
---|
76 | switches_(1), |
---|
77 | maximumTries_(-1) |
---|
78 | { |
---|
79 | if (howOften < -1900) { |
---|
80 | setGlobalCuts(true); |
---|
81 | howOften += 2000; |
---|
82 | } else if (howOften < -900) { |
---|
83 | setGlobalCutsAtRoot(true); |
---|
84 | howOften += 1000; |
---|
85 | } |
---|
86 | model_ = model; |
---|
87 | generator_ = generator->clone(); |
---|
88 | generator_->refreshSolver(model_->solver()); |
---|
89 | setNeedsOptimalBasis(generator_->needsOptimalBasis()); |
---|
90 | whenCutGenerator_ = howOften; |
---|
91 | whenCutGeneratorInSub_ = howOftenInSub; |
---|
92 | switchOffIfLessThan_ = switchOffIfLessThan; |
---|
93 | if (name) |
---|
94 | generatorName_ = CoinStrdup(name); |
---|
95 | else |
---|
96 | generatorName_ = CoinStrdup("Unknown"); |
---|
97 | setNormal(normal); |
---|
98 | setAtSolution(atSolution); |
---|
99 | setWhenInfeasible(infeasible); |
---|
100 | } |
---|
101 | |
---|
102 | // Copy constructor |
---|
103 | CbcCutGenerator::CbcCutGenerator ( const CbcCutGenerator & rhs) |
---|
104 | { |
---|
105 | model_ = rhs.model_; |
---|
106 | generator_ = rhs.generator_->clone(); |
---|
107 | //generator_->refreshSolver(model_->solver()); |
---|
108 | whenCutGenerator_ = rhs.whenCutGenerator_; |
---|
109 | whenCutGeneratorInSub_ = rhs.whenCutGeneratorInSub_; |
---|
110 | switchOffIfLessThan_ = rhs.switchOffIfLessThan_; |
---|
111 | depthCutGenerator_ = rhs.depthCutGenerator_; |
---|
112 | depthCutGeneratorInSub_ = rhs.depthCutGeneratorInSub_; |
---|
113 | generatorName_ = CoinStrdup(rhs.generatorName_); |
---|
114 | switches_ = rhs.switches_; |
---|
115 | maximumTries_ = rhs.maximumTries_; |
---|
116 | timeInCutGenerator_ = rhs.timeInCutGenerator_; |
---|
117 | savedCuts_ = rhs.savedCuts_; |
---|
118 | inaccuracy_ = rhs.inaccuracy_; |
---|
119 | numberTimes_ = rhs.numberTimes_; |
---|
120 | numberCuts_ = rhs.numberCuts_; |
---|
121 | numberElements_ = rhs.numberElements_; |
---|
122 | numberColumnCuts_ = rhs.numberColumnCuts_; |
---|
123 | numberCutsActive_ = rhs.numberCutsActive_; |
---|
124 | numberCutsAtRoot_ = rhs.numberCutsAtRoot_; |
---|
125 | numberActiveCutsAtRoot_ = rhs.numberActiveCutsAtRoot_; |
---|
126 | numberShortCutsAtRoot_ = rhs.numberShortCutsAtRoot_; |
---|
127 | } |
---|
128 | |
---|
129 | // Assignment operator |
---|
130 | CbcCutGenerator & |
---|
131 | CbcCutGenerator::operator=( const CbcCutGenerator & rhs) |
---|
132 | { |
---|
133 | if (this != &rhs) { |
---|
134 | delete generator_; |
---|
135 | free(generatorName_); |
---|
136 | model_ = rhs.model_; |
---|
137 | generator_ = rhs.generator_->clone(); |
---|
138 | generator_->refreshSolver(model_->solver()); |
---|
139 | whenCutGenerator_ = rhs.whenCutGenerator_; |
---|
140 | whenCutGeneratorInSub_ = rhs.whenCutGeneratorInSub_; |
---|
141 | switchOffIfLessThan_ = rhs.switchOffIfLessThan_; |
---|
142 | depthCutGenerator_ = rhs.depthCutGenerator_; |
---|
143 | depthCutGeneratorInSub_ = rhs.depthCutGeneratorInSub_; |
---|
144 | generatorName_ = CoinStrdup(rhs.generatorName_); |
---|
145 | switches_ = rhs.switches_; |
---|
146 | maximumTries_ = rhs.maximumTries_; |
---|
147 | timeInCutGenerator_ = rhs.timeInCutGenerator_; |
---|
148 | savedCuts_ = rhs.savedCuts_; |
---|
149 | inaccuracy_ = rhs.inaccuracy_; |
---|
150 | numberTimes_ = rhs.numberTimes_; |
---|
151 | numberCuts_ = rhs.numberCuts_; |
---|
152 | numberElements_ = rhs.numberElements_; |
---|
153 | numberColumnCuts_ = rhs.numberColumnCuts_; |
---|
154 | numberCutsActive_ = rhs.numberCutsActive_; |
---|
155 | numberCutsAtRoot_ = rhs.numberCutsAtRoot_; |
---|
156 | numberActiveCutsAtRoot_ = rhs.numberActiveCutsAtRoot_; |
---|
157 | numberShortCutsAtRoot_ = rhs.numberShortCutsAtRoot_; |
---|
158 | } |
---|
159 | return *this; |
---|
160 | } |
---|
161 | |
---|
162 | // Destructor |
---|
163 | CbcCutGenerator::~CbcCutGenerator () |
---|
164 | { |
---|
165 | free(generatorName_); |
---|
166 | delete generator_; |
---|
167 | } |
---|
168 | |
---|
169 | /* This is used to refresh any inforamtion. |
---|
170 | It also refreshes the solver in the cut generator |
---|
171 | in case generator wants to do some work |
---|
172 | */ |
---|
173 | void |
---|
174 | CbcCutGenerator::refreshModel(CbcModel * model) |
---|
175 | { |
---|
176 | model_ = model; |
---|
177 | // added test - helps if generator not thread safe |
---|
178 | if (whenCutGenerator_!=-100) |
---|
179 | generator_->refreshSolver(model_->solver()); |
---|
180 | } |
---|
181 | /* Generate cuts for the model data contained in si. |
---|
182 | The generated cuts are inserted into and returned in the |
---|
183 | collection of cuts cs. |
---|
184 | */ |
---|
185 | bool |
---|
186 | CbcCutGenerator::generateCuts( OsiCuts & cs , int fullScan, OsiSolverInterface * solver, CbcNode * node) |
---|
187 | { |
---|
188 | /* |
---|
189 | Make some decisions about whether we'll generate cuts. First convert |
---|
190 | whenCutGenerator_ to a set of canonical values for comparison to the node |
---|
191 | count. |
---|
192 | |
---|
193 | 0 < mod 1000000, with a result of 0 forced to 1 |
---|
194 | -99 <= <= 0 convert to 1 |
---|
195 | -100 = Off, period |
---|
196 | */ |
---|
197 | int depth; |
---|
198 | if (node) |
---|
199 | depth = node->depth(); |
---|
200 | else |
---|
201 | depth = 0; |
---|
202 | int howOften = whenCutGenerator_; |
---|
203 | if (dynamic_cast<CglProbing*>(generator_)) { |
---|
204 | if (howOften == -100 && model_->doCutsNow(3)) { |
---|
205 | howOften = 1; // do anyway |
---|
206 | } |
---|
207 | } |
---|
208 | if (howOften == -100) |
---|
209 | return false; |
---|
210 | int pass = model_->getCurrentPassNumber() - 1; |
---|
211 | if (maximumTries_>0) { |
---|
212 | // howOften means what it says |
---|
213 | if ((pass%howOften)!=0||depth) |
---|
214 | return false; |
---|
215 | else |
---|
216 | howOften=1; |
---|
217 | } |
---|
218 | if (howOften > 0) |
---|
219 | howOften = howOften % 1000000; |
---|
220 | else |
---|
221 | howOften = 1; |
---|
222 | if (!howOften) |
---|
223 | howOften = 1; |
---|
224 | bool returnCode = false; |
---|
225 | //OsiSolverInterface * solver = model_->solver(); |
---|
226 | // Reset cuts on first pass |
---|
227 | if (!pass) |
---|
228 | savedCuts_ = OsiCuts(); |
---|
229 | /* |
---|
230 | Determine if we should generate cuts based on node count. |
---|
231 | */ |
---|
232 | bool doThis = (model_->getNodeCount() % howOften) == 0; |
---|
233 | /* |
---|
234 | If the user has provided a depth specification, it will override the node |
---|
235 | count specification. |
---|
236 | */ |
---|
237 | if (depthCutGenerator_ > 0) { |
---|
238 | doThis = (depth % depthCutGenerator_) == 0; |
---|
239 | if (depth < depthCutGenerator_) |
---|
240 | doThis = true; // and also at top of tree |
---|
241 | } |
---|
242 | /* |
---|
243 | A few magic numbers ... |
---|
244 | |
---|
245 | The distinction between -100 and 100 for howOften is that we can override 100 |
---|
246 | with fullScan. -100 means no cuts, period. As does the magic number -200 for |
---|
247 | whenCutGeneratorInSub_. |
---|
248 | */ |
---|
249 | |
---|
250 | // But turn off if 100 |
---|
251 | if (howOften == 100) |
---|
252 | doThis = false; |
---|
253 | // Switch off if special setting |
---|
254 | if (whenCutGeneratorInSub_ == -200 && model_->parentModel()) { |
---|
255 | fullScan = 0; |
---|
256 | doThis = false; |
---|
257 | } |
---|
258 | if (fullScan || doThis) { |
---|
259 | CoinThreadRandom * randomNumberGenerator = NULL; |
---|
260 | #ifdef COIN_HAS_CLP |
---|
261 | { |
---|
262 | OsiClpSolverInterface * clpSolver |
---|
263 | = dynamic_cast<OsiClpSolverInterface *> (solver); |
---|
264 | if (clpSolver) |
---|
265 | randomNumberGenerator = |
---|
266 | clpSolver->getModelPtr()->randomNumberGenerator(); |
---|
267 | } |
---|
268 | #endif |
---|
269 | double time1 = 0.0; |
---|
270 | if (timing()) |
---|
271 | time1 = CoinCpuTime(); |
---|
272 | //#define CBC_DEBUG |
---|
273 | int numberRowCutsBefore = cs.sizeRowCuts() ; |
---|
274 | int numberColumnCutsBefore = cs.sizeColCuts() ; |
---|
275 | #ifdef JJF_ZERO |
---|
276 | int cutsBefore = cs.sizeCuts(); |
---|
277 | #endif |
---|
278 | CglTreeInfo info; |
---|
279 | info.level = depth; |
---|
280 | info.pass = pass; |
---|
281 | info.formulation_rows = model_->numberRowsAtContinuous(); |
---|
282 | info.inTree = node != NULL; |
---|
283 | info.randomNumberGenerator = randomNumberGenerator; |
---|
284 | info.options = (globalCutsAtRoot()) ? 8 : 0; |
---|
285 | if (ineffectualCuts()) |
---|
286 | info.options |= 32; |
---|
287 | if (globalCuts()) |
---|
288 | info.options |= 16; |
---|
289 | if (fullScan < 0) |
---|
290 | info.options |= 128; |
---|
291 | if (whetherInMustCallAgainMode()) |
---|
292 | info.options |= 1024; |
---|
293 | // See if we want alternate set of cuts |
---|
294 | if ((model_->moreSpecialOptions()&16384) != 0) |
---|
295 | info.options |= 256; |
---|
296 | if (model_->parentModel()) |
---|
297 | info.options |= 512; |
---|
298 | // above had &&!model_->parentModel()&&depth<2) |
---|
299 | incrementNumberTimesEntered(); |
---|
300 | CglProbing* generator = |
---|
301 | dynamic_cast<CglProbing*>(generator_); |
---|
302 | //if (!depth&&!pass) |
---|
303 | //printf("Cut generator %s when %d\n",generatorName_,whenCutGenerator_); |
---|
304 | if (!generator) { |
---|
305 | // Pass across model information in case it could be useful |
---|
306 | //void * saveData = solver->getApplicationData(); |
---|
307 | //solver->setApplicationData(model_); |
---|
308 | generator_->generateCuts(*solver, cs, info); |
---|
309 | //solver->setApplicationData(saveData); |
---|
310 | } else { |
---|
311 | // Probing - return tight column bounds |
---|
312 | CglTreeProbingInfo * info2 = model_->probingInfo(); |
---|
313 | bool doCuts = false; |
---|
314 | if (info2 && !depth) { |
---|
315 | info2->options = (globalCutsAtRoot()) ? 8 : 0; |
---|
316 | info2->level = depth; |
---|
317 | info2->pass = pass; |
---|
318 | info2->formulation_rows = model_->numberRowsAtContinuous(); |
---|
319 | info2->inTree = node != NULL; |
---|
320 | info2->randomNumberGenerator = randomNumberGenerator; |
---|
321 | generator->generateCutsAndModify(*solver, cs, info2); |
---|
322 | doCuts = true; |
---|
323 | } else if (depth) { |
---|
324 | /* The idea behind this is that probing may work in a different |
---|
325 | way deep in tree. So every now and then try various |
---|
326 | combinations to see what works. |
---|
327 | */ |
---|
328 | #define TRY_NOW_AND_THEN |
---|
329 | #ifdef TRY_NOW_AND_THEN |
---|
330 | if ((numberTimes_ == 200 || (numberTimes_ > 200 && (numberTimes_ % 2000) == 0)) |
---|
331 | && !model_->parentModel() && info.formulation_rows > 200) { |
---|
332 | /* In tree, every now and then try various combinations |
---|
333 | maxStack, maxProbe (last 5 digits) |
---|
334 | 123 is special and means CglProbing will try and |
---|
335 | be intelligent. |
---|
336 | */ |
---|
337 | int test[] = { |
---|
338 | 100123, |
---|
339 | 199999, |
---|
340 | 200123, |
---|
341 | 299999, |
---|
342 | 500123, |
---|
343 | 599999, |
---|
344 | 1000123, |
---|
345 | 1099999, |
---|
346 | 2000123, |
---|
347 | 2099999 |
---|
348 | }; |
---|
349 | int n = static_cast<int> (sizeof(test) / sizeof(int)); |
---|
350 | int saveStack = generator->getMaxLook(); |
---|
351 | int saveNumber = generator->getMaxProbe(); |
---|
352 | int kr1 = 0; |
---|
353 | int kc1 = 0; |
---|
354 | int bestStackTree = -1; |
---|
355 | int bestNumberTree = -1; |
---|
356 | for (int i = 0; i < n; i++) { |
---|
357 | //OsiCuts cs2 = cs; |
---|
358 | int stack = test[i] / 100000; |
---|
359 | int number = test[i] - 100000 * stack; |
---|
360 | generator->setMaxLook(stack); |
---|
361 | generator->setMaxProbe(number); |
---|
362 | int numberRowCutsBefore = cs.sizeRowCuts() ; |
---|
363 | int numberColumnCutsBefore = cs.sizeColCuts() ; |
---|
364 | generator_->generateCuts(*solver, cs, info); |
---|
365 | int numberRowCuts = cs.sizeRowCuts() - numberRowCutsBefore ; |
---|
366 | int numberColumnCuts = cs.sizeColCuts() - numberColumnCutsBefore ; |
---|
367 | #ifdef CLP_INVESTIGATE |
---|
368 | if (numberRowCuts < kr1 || numberColumnCuts < kc1) |
---|
369 | printf("Odd "); |
---|
370 | #endif |
---|
371 | if (numberRowCuts > kr1 || numberColumnCuts > kc1) { |
---|
372 | #ifdef CLP_INVESTIGATE |
---|
373 | printf("*** "); |
---|
374 | #endif |
---|
375 | kr1 = numberRowCuts; |
---|
376 | kc1 = numberColumnCuts; |
---|
377 | bestStackTree = stack; |
---|
378 | bestNumberTree = number; |
---|
379 | doCuts = true; |
---|
380 | } |
---|
381 | #ifdef CLP_INVESTIGATE |
---|
382 | printf("maxStack %d number %d gives %d row cuts and %d column cuts\n", |
---|
383 | stack, number, numberRowCuts, numberColumnCuts); |
---|
384 | #endif |
---|
385 | } |
---|
386 | generator->setMaxLook(saveStack); |
---|
387 | generator->setMaxProbe(saveNumber); |
---|
388 | if (bestStackTree > 0) { |
---|
389 | generator->setMaxLook(bestStackTree); |
---|
390 | generator->setMaxProbe(bestNumberTree); |
---|
391 | #ifdef CLP_INVESTIGATE |
---|
392 | printf("RRNumber %d -> %d, stack %d -> %d\n", |
---|
393 | saveNumber, bestNumberTree, saveStack, bestStackTree); |
---|
394 | #endif |
---|
395 | } else { |
---|
396 | // no good |
---|
397 | generator->setMaxLook(0); |
---|
398 | #ifdef CLP_INVESTIGATE |
---|
399 | printf("RRSwitching off number %d -> %d, stack %d -> %d\n", |
---|
400 | saveNumber, saveNumber, saveStack, 1); |
---|
401 | #endif |
---|
402 | } |
---|
403 | } |
---|
404 | #endif |
---|
405 | if (generator->getMaxLook() > 0 && !doCuts) { |
---|
406 | generator->generateCutsAndModify(*solver, cs, &info); |
---|
407 | doCuts = true; |
---|
408 | } |
---|
409 | } else { |
---|
410 | // at root - don't always do |
---|
411 | if (pass < 15 || (pass&1) == 0) { |
---|
412 | generator->generateCutsAndModify(*solver, cs, &info); |
---|
413 | doCuts = true; |
---|
414 | } |
---|
415 | } |
---|
416 | if (doCuts && generator->tightLower()) { |
---|
417 | // probing may have tightened bounds - check |
---|
418 | const double * tightLower = generator->tightLower(); |
---|
419 | const double * lower = solver->getColLower(); |
---|
420 | const double * tightUpper = generator->tightUpper(); |
---|
421 | const double * upper = solver->getColUpper(); |
---|
422 | const double * solution = solver->getColSolution(); |
---|
423 | int j; |
---|
424 | int numberColumns = solver->getNumCols(); |
---|
425 | double primalTolerance = 1.0e-8; |
---|
426 | const char * tightenBounds = generator->tightenBounds(); |
---|
427 | #ifdef CGL_DEBUG |
---|
428 | const OsiRowCutDebugger * debugger = solver->getRowCutDebugger(); |
---|
429 | if (debugger && debugger->onOptimalPath(*solver)) { |
---|
430 | printf("On optimal path CbcCut\n"); |
---|
431 | int nCols = solver->getNumCols(); |
---|
432 | int i; |
---|
433 | const double * optimal = debugger->optimalSolution(); |
---|
434 | const double * objective = solver->getObjCoefficients(); |
---|
435 | double objval1 = 0.0, objval2 = 0.0; |
---|
436 | for (i = 0; i < nCols; i++) { |
---|
437 | #if CGL_DEBUG>1 |
---|
438 | printf("%d %g %g %g %g\n", i, lower[i], solution[i], upper[i], optimal[i]); |
---|
439 | #endif |
---|
440 | objval1 += solution[i] * objective[i]; |
---|
441 | objval2 += optimal[i] * objective[i]; |
---|
442 | assert(optimal[i] >= lower[i] - 1.0e-5 && optimal[i] <= upper[i] + 1.0e-5); |
---|
443 | assert(optimal[i] >= tightLower[i] - 1.0e-5 && optimal[i] <= tightUpper[i] + 1.0e-5); |
---|
444 | } |
---|
445 | printf("current obj %g, integer %g\n", objval1, objval2); |
---|
446 | } |
---|
447 | #endif |
---|
448 | bool feasible = true; |
---|
449 | if ((model_->getThreadMode()&2) == 0) { |
---|
450 | for (j = 0; j < numberColumns; j++) { |
---|
451 | if (solver->isInteger(j)) { |
---|
452 | if (tightUpper[j] < upper[j]) { |
---|
453 | double nearest = floor(tightUpper[j] + 0.5); |
---|
454 | //assert (fabs(tightUpper[j]-nearest)<1.0e-5); may be infeasible |
---|
455 | solver->setColUpper(j, nearest); |
---|
456 | if (nearest < solution[j] - primalTolerance) |
---|
457 | returnCode = true; |
---|
458 | } |
---|
459 | if (tightLower[j] > lower[j]) { |
---|
460 | double nearest = floor(tightLower[j] + 0.5); |
---|
461 | //assert (fabs(tightLower[j]-nearest)<1.0e-5); may be infeasible |
---|
462 | solver->setColLower(j, nearest); |
---|
463 | if (nearest > solution[j] + primalTolerance) |
---|
464 | returnCode = true; |
---|
465 | } |
---|
466 | } else { |
---|
467 | if (upper[j] > lower[j]) { |
---|
468 | if (tightUpper[j] == tightLower[j]) { |
---|
469 | // fix |
---|
470 | //if (tightLower[j]!=lower[j]) |
---|
471 | solver->setColLower(j, tightLower[j]); |
---|
472 | //if (tightUpper[j]!=upper[j]) |
---|
473 | solver->setColUpper(j, tightUpper[j]); |
---|
474 | if (tightLower[j] > solution[j] + primalTolerance || |
---|
475 | tightUpper[j] < solution[j] - primalTolerance) |
---|
476 | returnCode = true; |
---|
477 | } else if (tightenBounds && tightenBounds[j]) { |
---|
478 | solver->setColLower(j, CoinMax(tightLower[j], lower[j])); |
---|
479 | solver->setColUpper(j, CoinMin(tightUpper[j], upper[j])); |
---|
480 | if (tightLower[j] > solution[j] + primalTolerance || |
---|
481 | tightUpper[j] < solution[j] - primalTolerance) |
---|
482 | returnCode = true; |
---|
483 | } |
---|
484 | } |
---|
485 | } |
---|
486 | if (upper[j] < lower[j] - 1.0e-3) { |
---|
487 | feasible = false; |
---|
488 | break; |
---|
489 | } |
---|
490 | } |
---|
491 | } else { |
---|
492 | CoinPackedVector lbs; |
---|
493 | CoinPackedVector ubs; |
---|
494 | int numberChanged = 0; |
---|
495 | bool ifCut = false; |
---|
496 | for (j = 0; j < numberColumns; j++) { |
---|
497 | if (solver->isInteger(j)) { |
---|
498 | if (tightUpper[j] < upper[j]) { |
---|
499 | double nearest = floor(tightUpper[j] + 0.5); |
---|
500 | //assert (fabs(tightUpper[j]-nearest)<1.0e-5); may be infeasible |
---|
501 | ubs.insert(j, nearest); |
---|
502 | numberChanged++; |
---|
503 | if (nearest < solution[j] - primalTolerance) |
---|
504 | ifCut = true; |
---|
505 | } |
---|
506 | if (tightLower[j] > lower[j]) { |
---|
507 | double nearest = floor(tightLower[j] + 0.5); |
---|
508 | //assert (fabs(tightLower[j]-nearest)<1.0e-5); may be infeasible |
---|
509 | lbs.insert(j, nearest); |
---|
510 | numberChanged++; |
---|
511 | if (nearest > solution[j] + primalTolerance) |
---|
512 | ifCut = true; |
---|
513 | } |
---|
514 | } else { |
---|
515 | if (upper[j] > lower[j]) { |
---|
516 | if (tightUpper[j] == tightLower[j]) { |
---|
517 | // fix |
---|
518 | lbs.insert(j, tightLower[j]); |
---|
519 | ubs.insert(j, tightUpper[j]); |
---|
520 | if (tightLower[j] > solution[j] + primalTolerance || |
---|
521 | tightUpper[j] < solution[j] - primalTolerance) |
---|
522 | ifCut = true; |
---|
523 | } else if (tightenBounds && tightenBounds[j]) { |
---|
524 | lbs.insert(j, CoinMax(tightLower[j], lower[j])); |
---|
525 | ubs.insert(j, CoinMin(tightUpper[j], upper[j])); |
---|
526 | if (tightLower[j] > solution[j] + primalTolerance || |
---|
527 | tightUpper[j] < solution[j] - primalTolerance) |
---|
528 | ifCut = true; |
---|
529 | } |
---|
530 | } |
---|
531 | } |
---|
532 | if (upper[j] < lower[j] - 1.0e-3) { |
---|
533 | feasible = false; |
---|
534 | break; |
---|
535 | } |
---|
536 | } |
---|
537 | if (numberChanged) { |
---|
538 | OsiColCut cc; |
---|
539 | cc.setUbs(ubs); |
---|
540 | cc.setLbs(lbs); |
---|
541 | if (ifCut) { |
---|
542 | cc.setEffectiveness(100.0); |
---|
543 | } else { |
---|
544 | cc.setEffectiveness(1.0e-5); |
---|
545 | } |
---|
546 | cs.insert(cc); |
---|
547 | } |
---|
548 | } |
---|
549 | if (!feasible) { |
---|
550 | // not feasible -add infeasible cut |
---|
551 | OsiRowCut rc; |
---|
552 | rc.setLb(COIN_DBL_MAX); |
---|
553 | rc.setUb(0.0); |
---|
554 | cs.insert(rc); |
---|
555 | } |
---|
556 | } |
---|
557 | //if (!solver->basisIsAvailable()) |
---|
558 | //returnCode=true; |
---|
559 | if (!returnCode) { |
---|
560 | // bounds changed but still optimal |
---|
561 | #ifdef COIN_HAS_CLP |
---|
562 | OsiClpSolverInterface * clpSolver |
---|
563 | = dynamic_cast<OsiClpSolverInterface *> (solver); |
---|
564 | if (clpSolver) { |
---|
565 | clpSolver->setLastAlgorithm(2); |
---|
566 | } |
---|
567 | #endif |
---|
568 | } |
---|
569 | #ifdef JJF_ZERO |
---|
570 | // Pass across info to pseudocosts |
---|
571 | char * mark = new char[numberColumns]; |
---|
572 | memset(mark, 0, numberColumns); |
---|
573 | int nLook = generator->numberThisTime(); |
---|
574 | const int * lookedAt = generator->lookedAt(); |
---|
575 | const int * fixedDown = generator->fixedDown(); |
---|
576 | const int * fixedUp = generator->fixedUp(); |
---|
577 | for (j = 0; j < nLook; j++) |
---|
578 | mark[lookedAt[j]] = 1; |
---|
579 | int numberObjects = model_->numberObjects(); |
---|
580 | for (int i = 0; i < numberObjects; i++) { |
---|
581 | CbcSimpleIntegerDynamicPseudoCost * obj1 = |
---|
582 | dynamic_cast <CbcSimpleIntegerDynamicPseudoCost *>(model_->modifiableObject(i)) ; |
---|
583 | if (obj1) { |
---|
584 | int iColumn = obj1->columnNumber(); |
---|
585 | if (mark[iColumn]) |
---|
586 | obj1->setProbingInformation(fixedDown[iColumn], fixedUp[iColumn]); |
---|
587 | } |
---|
588 | } |
---|
589 | delete [] mark; |
---|
590 | #endif |
---|
591 | } |
---|
592 | CbcCutModifier * modifier = model_->cutModifier(); |
---|
593 | if (modifier) { |
---|
594 | int numberRowCutsAfter = cs.sizeRowCuts() ; |
---|
595 | int k ; |
---|
596 | int nOdd = 0; |
---|
597 | //const OsiSolverInterface * solver = model_->solver(); |
---|
598 | for (k = numberRowCutsAfter - 1; k >= numberRowCutsBefore; k--) { |
---|
599 | OsiRowCut & thisCut = cs.rowCut(k) ; |
---|
600 | int returnCode = modifier->modify(solver, thisCut); |
---|
601 | if (returnCode) { |
---|
602 | nOdd++; |
---|
603 | if (returnCode == 3) |
---|
604 | cs.eraseRowCut(k); |
---|
605 | } |
---|
606 | } |
---|
607 | if (nOdd) |
---|
608 | COIN_DETAIL_PRINT(printf("Cut generator %s produced %d cuts of which %d were modified\n", |
---|
609 | generatorName_, numberRowCutsAfter - numberRowCutsBefore, nOdd)); |
---|
610 | } |
---|
611 | { |
---|
612 | // make all row cuts without test for duplicate |
---|
613 | int numberRowCutsAfter = cs.sizeRowCuts() ; |
---|
614 | int k ; |
---|
615 | #ifdef CGL_DEBUG |
---|
616 | const OsiRowCutDebugger * debugger = solver->getRowCutDebugger(); |
---|
617 | #endif |
---|
618 | for (k = numberRowCutsBefore; k < numberRowCutsAfter; k++) { |
---|
619 | OsiRowCut * thisCut = cs.rowCutPtr(k) ; |
---|
620 | #ifdef CGL_DEBUG |
---|
621 | if (debugger && debugger->onOptimalPath(*solver)) { |
---|
622 | assert(!debugger->invalidCut(*thisCut)); |
---|
623 | if(debugger->invalidCut(*thisCut)) |
---|
624 | abort(); |
---|
625 | } |
---|
626 | #endif |
---|
627 | thisCut->mutableRow().setTestForDuplicateIndex(false); |
---|
628 | } |
---|
629 | } |
---|
630 | // Add in saved cuts if violated |
---|
631 | if (false && !depth) { |
---|
632 | const double * solution = solver->getColSolution(); |
---|
633 | double primalTolerance = 1.0e-7; |
---|
634 | int numberCuts = savedCuts_.sizeRowCuts() ; |
---|
635 | for (int k = numberCuts - 1; k >= 0; k--) { |
---|
636 | const OsiRowCut * thisCut = savedCuts_.rowCutPtr(k) ; |
---|
637 | double sum = 0.0; |
---|
638 | int n = thisCut->row().getNumElements(); |
---|
639 | const int * column = thisCut->row().getIndices(); |
---|
640 | const double * element = thisCut->row().getElements(); |
---|
641 | assert (n); |
---|
642 | for (int i = 0; i < n; i++) { |
---|
643 | double value = element[i]; |
---|
644 | sum += value * solution[column[i]]; |
---|
645 | } |
---|
646 | if (sum > thisCut->ub() + primalTolerance) { |
---|
647 | sum = sum - thisCut->ub(); |
---|
648 | } else if (sum < thisCut->lb() - primalTolerance) { |
---|
649 | sum = thisCut->lb() - sum; |
---|
650 | } else { |
---|
651 | sum = 0.0; |
---|
652 | } |
---|
653 | if (sum) { |
---|
654 | // add to candidates and take out here |
---|
655 | cs.insert(*thisCut); |
---|
656 | savedCuts_.eraseRowCut(k); |
---|
657 | } |
---|
658 | } |
---|
659 | } |
---|
660 | if (!atSolution()) { |
---|
661 | int numberRowCutsAfter = cs.sizeRowCuts() ; |
---|
662 | int k ; |
---|
663 | int nEls = 0; |
---|
664 | int nCuts = numberRowCutsAfter - numberRowCutsBefore; |
---|
665 | // Remove NULL cuts! |
---|
666 | int nNull = 0; |
---|
667 | const double * solution = solver->getColSolution(); |
---|
668 | bool feasible = true; |
---|
669 | double primalTolerance = 1.0e-7; |
---|
670 | int shortCut = (depth) ? -1 : generator_->maximumLengthOfCutInTree(); |
---|
671 | for (k = numberRowCutsAfter - 1; k >= numberRowCutsBefore; k--) { |
---|
672 | const OsiRowCut * thisCut = cs.rowCutPtr(k) ; |
---|
673 | double sum = 0.0; |
---|
674 | if (thisCut->lb() <= thisCut->ub()) { |
---|
675 | int n = thisCut->row().getNumElements(); |
---|
676 | if (n <= shortCut) |
---|
677 | numberShortCutsAtRoot_++; |
---|
678 | const int * column = thisCut->row().getIndices(); |
---|
679 | const double * element = thisCut->row().getElements(); |
---|
680 | if (n <= 0) { |
---|
681 | // infeasible cut - give up |
---|
682 | feasible = false; |
---|
683 | break; |
---|
684 | } |
---|
685 | nEls += n; |
---|
686 | for (int i = 0; i < n; i++) { |
---|
687 | double value = element[i]; |
---|
688 | sum += value * solution[column[i]]; |
---|
689 | } |
---|
690 | if (sum > thisCut->ub() + primalTolerance) { |
---|
691 | sum = sum - thisCut->ub(); |
---|
692 | } else if (sum < thisCut->lb() - primalTolerance) { |
---|
693 | sum = thisCut->lb() - sum; |
---|
694 | } else { |
---|
695 | sum = 0.0; |
---|
696 | cs.eraseRowCut(k); |
---|
697 | nNull++; |
---|
698 | } |
---|
699 | } |
---|
700 | } |
---|
701 | //if (nNull) |
---|
702 | //printf("%s has %d cuts and %d elements - %d null!\n",generatorName_, |
---|
703 | // nCuts,nEls,nNull); |
---|
704 | numberRowCutsAfter = cs.sizeRowCuts() ; |
---|
705 | nCuts = numberRowCutsAfter - numberRowCutsBefore; |
---|
706 | nEls = 0; |
---|
707 | for (k = numberRowCutsBefore; k < numberRowCutsAfter; k++) { |
---|
708 | const OsiRowCut * thisCut = cs.rowCutPtr(k) ; |
---|
709 | int n = thisCut->row().getNumElements(); |
---|
710 | nEls += n; |
---|
711 | } |
---|
712 | //printf("%s has %d cuts and %d elements\n",generatorName_, |
---|
713 | // nCuts,nEls); |
---|
714 | int nElsNow = solver->getMatrixByCol()->getNumElements(); |
---|
715 | int numberColumns = solver->getNumCols(); |
---|
716 | int numberRows = solver->getNumRows(); |
---|
717 | //double averagePerRow = static_cast<double>(nElsNow)/ |
---|
718 | //static_cast<double>(numberRows); |
---|
719 | int nAdd; |
---|
720 | int nAdd2; |
---|
721 | int nReasonable; |
---|
722 | if (!model_->parentModel() && depth < 2) { |
---|
723 | if (inaccuracy_ < 3) { |
---|
724 | nAdd = 10000; |
---|
725 | if (pass > 0 && numberColumns > -500) |
---|
726 | nAdd = CoinMin(nAdd, nElsNow + 2 * numberRows); |
---|
727 | } else { |
---|
728 | nAdd = 10000; |
---|
729 | if (pass > 0) |
---|
730 | nAdd = CoinMin(nAdd, nElsNow + 2 * numberRows); |
---|
731 | } |
---|
732 | nAdd2 = 5 * numberColumns; |
---|
733 | nReasonable = CoinMax(nAdd2, nElsNow / 8 + nAdd); |
---|
734 | if (!depth && !pass) { |
---|
735 | // allow more |
---|
736 | nAdd += nElsNow / 2; |
---|
737 | nAdd2 += nElsNow / 2; |
---|
738 | nReasonable += nElsNow / 2; |
---|
739 | } |
---|
740 | //if (!depth&&ineffectualCuts()) |
---|
741 | //nReasonable *= 2; |
---|
742 | } else { |
---|
743 | nAdd = 200; |
---|
744 | nAdd2 = 2 * numberColumns; |
---|
745 | nReasonable = CoinMax(nAdd2, nElsNow / 8 + nAdd); |
---|
746 | } |
---|
747 | //#define UNS_WEIGHT 0.1 |
---|
748 | #ifdef UNS_WEIGHT |
---|
749 | const double * colLower = solver->getColLower(); |
---|
750 | const double * colUpper = solver->getColUpper(); |
---|
751 | #endif |
---|
752 | if (/*nEls>CoinMax(nAdd2,nElsNow/8+nAdd)*/nCuts && feasible) { |
---|
753 | //printf("need to remove cuts\n"); |
---|
754 | // just add most effective |
---|
755 | #ifndef JJF_ONE |
---|
756 | int nDelete = nEls - nReasonable; |
---|
757 | |
---|
758 | nElsNow = nEls; |
---|
759 | double * sort = new double [nCuts]; |
---|
760 | int * which = new int [nCuts]; |
---|
761 | // For parallel cuts |
---|
762 | double * element2 = new double [numberColumns]; |
---|
763 | //#define USE_OBJECTIVE 2 |
---|
764 | #ifdef USE_OBJECTIVE |
---|
765 | const double *objective = solver->getObjCoefficients() ; |
---|
766 | #if USE_OBJECTIVE>1 |
---|
767 | double objNorm = 0.0; |
---|
768 | for (int i = 0; i < numberColumns; i++) |
---|
769 | objNorm += objective[i] * objective[i]; |
---|
770 | if (objNorm) |
---|
771 | objNorm = 1.0 / sqrt(objNorm); |
---|
772 | else |
---|
773 | objNorm = 1.0; |
---|
774 | objNorm *= 0.01; // downgrade |
---|
775 | #endif |
---|
776 | #endif |
---|
777 | CoinZeroN(element2, numberColumns); |
---|
778 | for (k = numberRowCutsBefore; k < numberRowCutsAfter; k++) { |
---|
779 | const OsiRowCut * thisCut = cs.rowCutPtr(k) ; |
---|
780 | double sum = 0.0; |
---|
781 | if (thisCut->lb() <= thisCut->ub()) { |
---|
782 | int n = thisCut->row().getNumElements(); |
---|
783 | const int * column = thisCut->row().getIndices(); |
---|
784 | const double * element = thisCut->row().getElements(); |
---|
785 | assert (n); |
---|
786 | #ifdef UNS_WEIGHT |
---|
787 | double normU = 0.0; |
---|
788 | double norm = 1.0e-3; |
---|
789 | int nU = 0; |
---|
790 | for (int i = 0; i < n; i++) { |
---|
791 | double value = element[i]; |
---|
792 | int iColumn = column[i]; |
---|
793 | double solValue = solution[iColumn]; |
---|
794 | sum += value * solValue; |
---|
795 | value *= value; |
---|
796 | norm += value; |
---|
797 | if (solValue > colLower[iColumn] + 1.0e-6 && |
---|
798 | solValue < colUpper[iColumn] - 1.0e-6) { |
---|
799 | normU += value; |
---|
800 | nU++; |
---|
801 | } |
---|
802 | } |
---|
803 | #ifdef JJF_ZERO |
---|
804 | int nS = n - nU; |
---|
805 | if (numberColumns > 20000) { |
---|
806 | if (nS > 50) { |
---|
807 | double ratio = 50.0 / nS; |
---|
808 | normU /= ratio; |
---|
809 | } |
---|
810 | } |
---|
811 | #endif |
---|
812 | norm += UNS_WEIGHT * (normU - norm); |
---|
813 | #else |
---|
814 | double norm = 1.0e-3; |
---|
815 | #ifdef USE_OBJECTIVE |
---|
816 | double obj = 0.0; |
---|
817 | #endif |
---|
818 | for (int i = 0; i < n; i++) { |
---|
819 | int iColumn = column[i]; |
---|
820 | double value = element[i]; |
---|
821 | sum += value * solution[iColumn]; |
---|
822 | norm += value * value; |
---|
823 | #ifdef USE_OBJECTIVE |
---|
824 | obj += value * objective[iColumn]; |
---|
825 | #endif |
---|
826 | } |
---|
827 | #endif |
---|
828 | if (sum > thisCut->ub()) { |
---|
829 | sum = sum - thisCut->ub(); |
---|
830 | } else if (sum < thisCut->lb()) { |
---|
831 | sum = thisCut->lb() - sum; |
---|
832 | } else { |
---|
833 | sum = 0.0; |
---|
834 | } |
---|
835 | #ifdef USE_OBJECTIVE |
---|
836 | if (sum) { |
---|
837 | #if USE_OBJECTIVE==1 |
---|
838 | obj = CoinMax(1.0e-6, fabs(obj)); |
---|
839 | norm = sqrt(obj * norm); |
---|
840 | //sum += fabs(obj)*invObjNorm; |
---|
841 | //printf("sum %g norm %g normobj %g invNorm %g mod %g\n", |
---|
842 | // sum,norm,obj,invObjNorm,obj*invObjNorm); |
---|
843 | // normalize |
---|
844 | sum /= sqrt(norm); |
---|
845 | #else |
---|
846 | // normalize |
---|
847 | norm = 1.0 / sqrt(norm); |
---|
848 | sum = (sum + objNorm * obj) * norm; |
---|
849 | #endif |
---|
850 | } |
---|
851 | #else |
---|
852 | // normalize |
---|
853 | sum /= sqrt(norm); |
---|
854 | #endif |
---|
855 | //sum /= pow(norm,0.3); |
---|
856 | // adjust for length |
---|
857 | //sum /= pow(reinterpret_cast<double>(n),0.2); |
---|
858 | //sum /= sqrt((double) n); |
---|
859 | // randomize |
---|
860 | //double randomNumber = |
---|
861 | //model_->randomNumberGenerator()->randomDouble(); |
---|
862 | //sum *= (0.5+randomNumber); |
---|
863 | } else { |
---|
864 | // keep |
---|
865 | sum = COIN_DBL_MAX; |
---|
866 | } |
---|
867 | sort[k-numberRowCutsBefore] = sum; |
---|
868 | which[k-numberRowCutsBefore] = k; |
---|
869 | } |
---|
870 | CoinSort_2(sort, sort + nCuts, which); |
---|
871 | // Now see which ones are too similar |
---|
872 | int nParallel = 0; |
---|
873 | double testValue = (depth > 1) ? 0.99 : 0.999999; |
---|
874 | for (k = 0; k < nCuts; k++) { |
---|
875 | int j = which[k]; |
---|
876 | const OsiRowCut * thisCut = cs.rowCutPtr(j) ; |
---|
877 | if (thisCut->lb() > thisCut->ub()) |
---|
878 | break; // cut is infeasible |
---|
879 | int n = thisCut->row().getNumElements(); |
---|
880 | const int * column = thisCut->row().getIndices(); |
---|
881 | const double * element = thisCut->row().getElements(); |
---|
882 | assert (n); |
---|
883 | double norm = 0.0; |
---|
884 | double lb = thisCut->lb(); |
---|
885 | double ub = thisCut->ub(); |
---|
886 | for (int i = 0; i < n; i++) { |
---|
887 | double value = element[i]; |
---|
888 | element2[column[i]] = value; |
---|
889 | norm += value * value; |
---|
890 | } |
---|
891 | int kkk = CoinMin(nCuts, k + 5); |
---|
892 | for (int kk = k + 1; kk < kkk; kk++) { |
---|
893 | int jj = which[kk]; |
---|
894 | const OsiRowCut * thisCut2 = cs.rowCutPtr(jj) ; |
---|
895 | if (thisCut2->lb() > thisCut2->ub()) |
---|
896 | break; // cut is infeasible |
---|
897 | int nB = thisCut2->row().getNumElements(); |
---|
898 | const int * columnB = thisCut2->row().getIndices(); |
---|
899 | const double * elementB = thisCut2->row().getElements(); |
---|
900 | assert (nB); |
---|
901 | double normB = 0.0; |
---|
902 | double product = 0.0; |
---|
903 | for (int i = 0; i < nB; i++) { |
---|
904 | double value = elementB[i]; |
---|
905 | normB += value * value; |
---|
906 | product += value * element2[columnB[i]]; |
---|
907 | } |
---|
908 | if (product > 0.0 && product*product > testValue*norm*normB) { |
---|
909 | bool parallel = true; |
---|
910 | double lbB = thisCut2->lb(); |
---|
911 | double ubB = thisCut2->ub(); |
---|
912 | if ((lb < -1.0e20 && lbB > -1.0e20) || |
---|
913 | (lbB < -1.0e20 && lb > -1.0e20)) |
---|
914 | parallel = false; |
---|
915 | double tolerance; |
---|
916 | tolerance = CoinMax(fabs(lb), fabs(lbB)) + 1.0e-6; |
---|
917 | if (fabs(lb - lbB) > tolerance) |
---|
918 | parallel = false; |
---|
919 | if ((ub > 1.0e20 && ubB < 1.0e20) || |
---|
920 | (ubB > 1.0e20 && ub < 1.0e20)) |
---|
921 | parallel = false; |
---|
922 | tolerance = CoinMax(fabs(ub), fabs(ubB)) + 1.0e-6; |
---|
923 | if (fabs(ub - ubB) > tolerance) |
---|
924 | parallel = false; |
---|
925 | if (parallel) { |
---|
926 | nParallel++; |
---|
927 | sort[k] = 0.0; |
---|
928 | break; |
---|
929 | } |
---|
930 | } |
---|
931 | } |
---|
932 | for (int i = 0; i < n; i++) { |
---|
933 | element2[column[i]] = 0.0; |
---|
934 | } |
---|
935 | } |
---|
936 | delete [] element2; |
---|
937 | CoinSort_2(sort, sort + nCuts, which); |
---|
938 | k = 0; |
---|
939 | while (nDelete > 0 || !sort[k]) { |
---|
940 | int iCut = which[k]; |
---|
941 | const OsiRowCut * thisCut = cs.rowCutPtr(iCut) ; |
---|
942 | int n = thisCut->row().getNumElements(); |
---|
943 | // may be best, just to save if short |
---|
944 | if (false && n && sort[k]) { |
---|
945 | // add to saved cuts |
---|
946 | savedCuts_.insert(*thisCut); |
---|
947 | } |
---|
948 | nDelete -= n; |
---|
949 | k++; |
---|
950 | if (k >= nCuts) |
---|
951 | break; |
---|
952 | } |
---|
953 | std::sort(which, which + k); |
---|
954 | k--; |
---|
955 | for (; k >= 0; k--) { |
---|
956 | cs.eraseRowCut(which[k]); |
---|
957 | } |
---|
958 | delete [] sort; |
---|
959 | delete [] which; |
---|
960 | numberRowCutsAfter = cs.sizeRowCuts() ; |
---|
961 | #else |
---|
962 | double * norm = new double [nCuts]; |
---|
963 | int * which = new int [2*nCuts]; |
---|
964 | double * score = new double [nCuts]; |
---|
965 | double * ortho = new double [nCuts]; |
---|
966 | int nIn = 0; |
---|
967 | int nOut = nCuts; |
---|
968 | // For parallel cuts |
---|
969 | double * element2 = new double [numberColumns]; |
---|
970 | const double *objective = solver->getObjCoefficients() ; |
---|
971 | double objNorm = 0.0; |
---|
972 | for (int i = 0; i < numberColumns; i++) |
---|
973 | objNorm += objective[i] * objective[i]; |
---|
974 | if (objNorm) |
---|
975 | objNorm = 1.0 / sqrt(objNorm); |
---|
976 | else |
---|
977 | objNorm = 1.0; |
---|
978 | objNorm *= 0.1; // weight of 0.1 |
---|
979 | CoinZeroN(element2, numberColumns); |
---|
980 | int numberRowCuts = numberRowCutsAfter - numberRowCutsBefore; |
---|
981 | int iBest = -1; |
---|
982 | double best = 0.0; |
---|
983 | int nPossible = 0; |
---|
984 | double testValue = (depth > 1) ? 0.7 : 0.5; |
---|
985 | for (k = 0; k < numberRowCuts; k++) { |
---|
986 | const OsiRowCut * thisCut = cs.rowCutPtr(k + numberRowCutsBefore) ; |
---|
987 | double sum = 0.0; |
---|
988 | if (thisCut->lb() <= thisCut->ub()) { |
---|
989 | int n = thisCut->row().getNumElements(); |
---|
990 | const int * column = thisCut->row().getIndices(); |
---|
991 | const double * element = thisCut->row().getElements(); |
---|
992 | assert (n); |
---|
993 | double normThis = 1.0e-6; |
---|
994 | double obj = 0.0; |
---|
995 | for (int i = 0; i < n; i++) { |
---|
996 | int iColumn = column[i]; |
---|
997 | double value = element[i]; |
---|
998 | sum += value * solution[iColumn]; |
---|
999 | normThis += value * value; |
---|
1000 | obj += value * objective[iColumn]; |
---|
1001 | } |
---|
1002 | if (sum > thisCut->ub()) { |
---|
1003 | sum = sum - thisCut->ub(); |
---|
1004 | } else if (sum < thisCut->lb()) { |
---|
1005 | sum = thisCut->lb() - sum; |
---|
1006 | } else { |
---|
1007 | sum = 0.0; |
---|
1008 | } |
---|
1009 | if (sum) { |
---|
1010 | normThis = 1.0 / sqrt(normThis); |
---|
1011 | norm[k] = normThis; |
---|
1012 | sum *= normThis; |
---|
1013 | obj *= normThis; |
---|
1014 | score[k] = sum + obj * objNorm; |
---|
1015 | ortho[k] = 1.0; |
---|
1016 | } |
---|
1017 | } else { |
---|
1018 | // keep and discard others |
---|
1019 | nIn = 1; |
---|
1020 | which[0] = k; |
---|
1021 | for (int j = 0; j < numberRowCuts; j++) { |
---|
1022 | if (j != k) |
---|
1023 | which[nOut++] = j; |
---|
1024 | } |
---|
1025 | iBest = -1; |
---|
1026 | break; |
---|
1027 | } |
---|
1028 | if (sum) { |
---|
1029 | if (score[k] > best) { |
---|
1030 | best = score[k]; |
---|
1031 | iBest = nPossible; |
---|
1032 | } |
---|
1033 | which[nPossible++] = k; |
---|
1034 | } else { |
---|
1035 | which[nOut++] = k; |
---|
1036 | } |
---|
1037 | } |
---|
1038 | while (iBest >= 0) { |
---|
1039 | int kBest = which[iBest]; |
---|
1040 | int j = which[nIn]; |
---|
1041 | which[iBest] = j; |
---|
1042 | which[nIn++] = kBest; |
---|
1043 | const OsiRowCut * thisCut = cs.rowCutPtr(kBest + numberRowCutsBefore) ; |
---|
1044 | int n = thisCut->row().getNumElements(); |
---|
1045 | nReasonable -= n; |
---|
1046 | if (nReasonable <= 0) { |
---|
1047 | for (k = nIn; k < nPossible; k++) |
---|
1048 | which[nOut++] = which[k]; |
---|
1049 | break; |
---|
1050 | } |
---|
1051 | // Now see which ones are too similar and choose next |
---|
1052 | iBest = -1; |
---|
1053 | best = 0.0; |
---|
1054 | int nOld = nPossible; |
---|
1055 | nPossible = nIn; |
---|
1056 | const int * column = thisCut->row().getIndices(); |
---|
1057 | const double * element = thisCut->row().getElements(); |
---|
1058 | assert (n); |
---|
1059 | double normNew = norm[kBest]; |
---|
1060 | for (int i = 0; i < n; i++) { |
---|
1061 | double value = element[i]; |
---|
1062 | element2[column[i]] = value; |
---|
1063 | } |
---|
1064 | for (int j = nIn; j < nOld; j++) { |
---|
1065 | k = which[j]; |
---|
1066 | const OsiRowCut * thisCut2 = cs.rowCutPtr(k + numberRowCutsBefore) ; |
---|
1067 | int nB = thisCut2->row().getNumElements(); |
---|
1068 | const int * columnB = thisCut2->row().getIndices(); |
---|
1069 | const double * elementB = thisCut2->row().getElements(); |
---|
1070 | assert (nB); |
---|
1071 | double normB = norm[k]; |
---|
1072 | double product = 0.0; |
---|
1073 | for (int i = 0; i < nB; i++) { |
---|
1074 | double value = elementB[i]; |
---|
1075 | product += value * element2[columnB[i]]; |
---|
1076 | } |
---|
1077 | double orthoScore = 1.0 - product * normNew * normB; |
---|
1078 | if (orthoScore >= testValue) { |
---|
1079 | ortho[k] = CoinMin(orthoScore, ortho[k]); |
---|
1080 | double test = score[k] + ortho[k]; |
---|
1081 | if (test > best) { |
---|
1082 | best = score[k]; |
---|
1083 | iBest = nPossible; |
---|
1084 | } |
---|
1085 | which[nPossible++] = k; |
---|
1086 | } else { |
---|
1087 | which[nOut++] = k; |
---|
1088 | } |
---|
1089 | } |
---|
1090 | for (int i = 0; i < n; i++) { |
---|
1091 | element2[column[i]] = 0.0; |
---|
1092 | } |
---|
1093 | } |
---|
1094 | delete [] score; |
---|
1095 | delete [] ortho; |
---|
1096 | std::sort(which + nCuts, which + nOut); |
---|
1097 | k = nOut - 1; |
---|
1098 | for (; k >= nCuts; k--) { |
---|
1099 | cs.eraseRowCut(which[k] + numberRowCutsBefore); |
---|
1100 | } |
---|
1101 | delete [] norm; |
---|
1102 | delete [] which; |
---|
1103 | numberRowCutsAfter = cs.sizeRowCuts() ; |
---|
1104 | #endif |
---|
1105 | } |
---|
1106 | } |
---|
1107 | #ifdef CBC_DEBUG |
---|
1108 | { |
---|
1109 | int numberRowCutsAfter = cs.sizeRowCuts() ; |
---|
1110 | int k ; |
---|
1111 | int nBad = 0; |
---|
1112 | for (k = numberRowCutsBefore; k < numberRowCutsAfter; k++) { |
---|
1113 | OsiRowCut thisCut = cs.rowCut(k) ; |
---|
1114 | if (thisCut.lb() > thisCut.ub() || |
---|
1115 | thisCut.lb() > 1.0e8 || |
---|
1116 | thisCut.ub() < -1.0e8) |
---|
1117 | printf("cut from %s has bounds %g and %g!\n", |
---|
1118 | generatorName_, thisCut.lb(), thisCut.ub()); |
---|
1119 | if (thisCut.lb() <= thisCut.ub()) { |
---|
1120 | /* check size of elements. |
---|
1121 | We can allow smaller but this helps debug generators as it |
---|
1122 | is unsafe to have small elements */ |
---|
1123 | int n = thisCut.row().getNumElements(); |
---|
1124 | const int * column = thisCut.row().getIndices(); |
---|
1125 | const double * element = thisCut.row().getElements(); |
---|
1126 | assert (n); |
---|
1127 | for (int i = 0; i < n; i++) { |
---|
1128 | double value = element[i]; |
---|
1129 | if (fabs(value) <= 1.0e-12 || fabs(value) >= 1.0e20) |
---|
1130 | nBad++; |
---|
1131 | } |
---|
1132 | } |
---|
1133 | if (nBad) |
---|
1134 | printf("Cut generator %s produced %d cuts of which %d had tiny or large elements\n", |
---|
1135 | generatorName_, numberRowCutsAfter - numberRowCutsBefore, nBad); |
---|
1136 | } |
---|
1137 | } |
---|
1138 | #endif |
---|
1139 | int numberRowCutsAfter = cs.sizeRowCuts() ; |
---|
1140 | int numberColumnCutsAfter = cs.sizeColCuts() ; |
---|
1141 | if (numberRowCutsBefore < numberRowCutsAfter) { |
---|
1142 | for (int k = numberRowCutsBefore; k < numberRowCutsAfter; k++) { |
---|
1143 | OsiRowCut thisCut = cs.rowCut(k) ; |
---|
1144 | int n = thisCut.row().getNumElements(); |
---|
1145 | numberElements_ += n; |
---|
1146 | } |
---|
1147 | #ifdef JJF_ZERO |
---|
1148 | printf("generator %s generated %d row cuts\n", |
---|
1149 | generatorName_, numberRowCutsAfter - numberRowCutsBefore); |
---|
1150 | #endif |
---|
1151 | numberCuts_ += numberRowCutsAfter - numberRowCutsBefore; |
---|
1152 | } |
---|
1153 | if (numberColumnCutsBefore < numberColumnCutsAfter) { |
---|
1154 | #ifdef JJF_ZERO |
---|
1155 | printf("generator %s generated %d column cuts\n", |
---|
1156 | generatorName_, numberColumnCutsAfter - numberColumnCutsBefore); |
---|
1157 | #endif |
---|
1158 | numberColumnCuts_ += numberColumnCutsAfter - numberColumnCutsBefore; |
---|
1159 | } |
---|
1160 | if (timing()) |
---|
1161 | timeInCutGenerator_ += CoinCpuTime() - time1; |
---|
1162 | // switch off if first time and no good |
---|
1163 | if (node == NULL && !pass ) { |
---|
1164 | if (numberRowCutsAfter - numberRowCutsBefore |
---|
1165 | < switchOffIfLessThan_ /*&& numberCuts_ < switchOffIfLessThan_*/) { |
---|
1166 | // switch off |
---|
1167 | maximumTries_ = 0; |
---|
1168 | whenCutGenerator_=-100; |
---|
1169 | //whenCutGenerator_ = -100; |
---|
1170 | //whenCutGeneratorInSub_ = -200; |
---|
1171 | } |
---|
1172 | } |
---|
1173 | if (maximumTries_>0) { |
---|
1174 | maximumTries_--; |
---|
1175 | if (!maximumTries_) |
---|
1176 | whenCutGenerator_=-100; |
---|
1177 | } |
---|
1178 | } |
---|
1179 | return returnCode; |
---|
1180 | } |
---|
1181 | void |
---|
1182 | CbcCutGenerator::setHowOften(int howOften) |
---|
1183 | { |
---|
1184 | |
---|
1185 | if (howOften >= 1000000) { |
---|
1186 | // leave Probing every SCANCUTS_PROBING |
---|
1187 | howOften = howOften % 1000000; |
---|
1188 | CglProbing* generator = |
---|
1189 | dynamic_cast<CglProbing*>(generator_); |
---|
1190 | |
---|
1191 | if (generator && howOften > SCANCUTS_PROBING) |
---|
1192 | howOften = SCANCUTS_PROBING + 1000000; |
---|
1193 | else |
---|
1194 | howOften += 1000000; |
---|
1195 | } |
---|
1196 | whenCutGenerator_ = howOften; |
---|
1197 | } |
---|
1198 | void |
---|
1199 | CbcCutGenerator::setWhatDepth(int value) |
---|
1200 | { |
---|
1201 | depthCutGenerator_ = value; |
---|
1202 | } |
---|
1203 | void |
---|
1204 | CbcCutGenerator::setWhatDepthInSub(int value) |
---|
1205 | { |
---|
1206 | depthCutGeneratorInSub_ = value; |
---|
1207 | } |
---|
1208 | // Add in statistics from other |
---|
1209 | void |
---|
1210 | CbcCutGenerator::addStatistics(const CbcCutGenerator * other) |
---|
1211 | { |
---|
1212 | // Time in cut generator |
---|
1213 | timeInCutGenerator_ += other->timeInCutGenerator_; |
---|
1214 | // Number times cut generator entered |
---|
1215 | numberTimes_ += other->numberTimes_; |
---|
1216 | // Total number of cuts added |
---|
1217 | numberCuts_ += other->numberCuts_; |
---|
1218 | // Total number of elements added |
---|
1219 | numberElements_ += other->numberElements_; |
---|
1220 | // Total number of column cuts added |
---|
1221 | numberColumnCuts_ += other->numberColumnCuts_; |
---|
1222 | // Total number of cuts active after (at end of n cut passes at each node) |
---|
1223 | numberCutsActive_ += other->numberCutsActive_; |
---|
1224 | // Number of cuts generated at root |
---|
1225 | numberCutsAtRoot_ += other->numberCutsAtRoot_; |
---|
1226 | // Number of cuts active at root |
---|
1227 | numberActiveCutsAtRoot_ += other->numberActiveCutsAtRoot_; |
---|
1228 | // Number of short cuts at root |
---|
1229 | numberShortCutsAtRoot_ += other->numberShortCutsAtRoot_; |
---|
1230 | } |
---|
1231 | // Scale back statistics by factor |
---|
1232 | void |
---|
1233 | CbcCutGenerator::scaleBackStatistics(int factor) |
---|
1234 | { |
---|
1235 | // leave time |
---|
1236 | // Number times cut generator entered |
---|
1237 | numberTimes_ = (numberTimes_+factor-1)/factor; |
---|
1238 | // Total number of cuts added |
---|
1239 | numberCuts_ = (numberCuts_+factor-1)/factor; |
---|
1240 | // Total number of elements added |
---|
1241 | numberElements_ = (numberElements_+factor-1)/factor; |
---|
1242 | // Total number of column cuts added |
---|
1243 | numberColumnCuts_ = (numberColumnCuts_+factor-1)/factor; |
---|
1244 | // Total number of cuts active after (at end of n cut passes at each node) |
---|
1245 | numberCutsActive_ = (numberCutsActive_+factor-1)/factor; |
---|
1246 | // Number of cuts generated at root |
---|
1247 | numberCutsAtRoot_ = (numberCutsAtRoot_+factor-1)/factor; |
---|
1248 | // Number of cuts active at root |
---|
1249 | numberActiveCutsAtRoot_ = (numberActiveCutsAtRoot_+factor-1)/factor; |
---|
1250 | // Number of short cuts at root |
---|
1251 | numberShortCutsAtRoot_ = (numberShortCutsAtRoot_+factor-1)/factor; |
---|
1252 | } |
---|
1253 | // Create C++ lines to get to current state |
---|
1254 | void |
---|
1255 | CbcCutGenerator::generateTuning( FILE * fp) |
---|
1256 | { |
---|
1257 | fprintf(fp, "// Cbc tuning for generator %s\n", generatorName_); |
---|
1258 | fprintf(fp, " generator->setHowOften(%d);\n", whenCutGenerator_); |
---|
1259 | fprintf(fp, " generator->setSwitchOffIfLessThan(%d);\n", switchOffIfLessThan_); |
---|
1260 | fprintf(fp, " generator->setWhatDepth(%d);\n", depthCutGenerator_); |
---|
1261 | fprintf(fp, " generator->setInaccuracy(%d);\n", inaccuracy_); |
---|
1262 | if (timing()) |
---|
1263 | fprintf(fp, " generator->setTiming(true);\n"); |
---|
1264 | if (normal()) |
---|
1265 | fprintf(fp, " generator->setNormal(true);\n"); |
---|
1266 | if (atSolution()) |
---|
1267 | fprintf(fp, " generator->setAtSolution(true);\n"); |
---|
1268 | if (whenInfeasible()) |
---|
1269 | fprintf(fp, " generator->setWhenInfeasible(true);\n"); |
---|
1270 | if (needsOptimalBasis()) |
---|
1271 | fprintf(fp, " generator->setNeedsOptimalBasis(true);\n"); |
---|
1272 | if (mustCallAgain()) |
---|
1273 | fprintf(fp, " generator->setMustCallAgain(true);\n"); |
---|
1274 | if (whetherToUse()) |
---|
1275 | fprintf(fp, " generator->setWhetherToUse(true);\n"); |
---|
1276 | } |
---|
1277 | |
---|
1278 | |
---|
1279 | |
---|