1 | /* $Id: CbcCutGenerator.cpp 2151 2015-03-04 15:40:35Z 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 | //#define WEAKEN_CUTS 1 |
---|
619 | #ifdef WEAKEN_CUTS |
---|
620 | const double * lower = solver->getColLower(); |
---|
621 | const double * upper = solver->getColUpper(); |
---|
622 | const double * solution = solver->getColSolution(); |
---|
623 | #endif |
---|
624 | for (k = numberRowCutsBefore; k < numberRowCutsAfter; k++) { |
---|
625 | OsiRowCut * thisCut = cs.rowCutPtr(k) ; |
---|
626 | #ifdef WEAKEN_CUTS |
---|
627 | // weaken cut if coefficients not integer |
---|
628 | |
---|
629 | double lb=thisCut->lb(); |
---|
630 | double ub=thisCut->ub(); |
---|
631 | if (lb<-1.0e100||ub>1.0e100) { |
---|
632 | // normal cut |
---|
633 | CoinPackedVector rpv = thisCut->row(); |
---|
634 | const int n = rpv.getNumElements(); |
---|
635 | const int * indices = rpv.getIndices(); |
---|
636 | const double * elements = rpv.getElements(); |
---|
637 | double bound=0.0; |
---|
638 | double sum=0.0; |
---|
639 | bool integral=true; |
---|
640 | int nInteger=0; |
---|
641 | for (int k=0; k<n; k++) { |
---|
642 | double value = fabs(elements[k]); |
---|
643 | int column=indices[k]; |
---|
644 | sum += value; |
---|
645 | if (value!=floor(value+0.5)) |
---|
646 | integral=false; |
---|
647 | if (solver->isInteger(column)) { |
---|
648 | nInteger++; |
---|
649 | double largerBound = CoinMax(fabs(lower[column]), |
---|
650 | fabs(upper[column])); |
---|
651 | double solutionBound=fabs(solution[column])+10.0; |
---|
652 | bound += CoinMin(largerBound,solutionBound); |
---|
653 | } |
---|
654 | } |
---|
655 | #if WEAKEN_CUTS ==1 |
---|
656 | // leave if all 0-1 |
---|
657 | if (nInteger==bound) |
---|
658 | integral=true; |
---|
659 | #endif |
---|
660 | if (!integral) { |
---|
661 | double weakenBy=1.0e-7*(bound+sum); |
---|
662 | #if WEAKEN_CUTS>2 |
---|
663 | weakenBy *= 10.0; |
---|
664 | #endif |
---|
665 | if (lb<-1.0e100) |
---|
666 | thisCut->setUb(ub+weakenBy); |
---|
667 | else |
---|
668 | thisCut->setLb(lb-weakenBy); |
---|
669 | } |
---|
670 | } |
---|
671 | #endif |
---|
672 | #ifdef CGL_DEBUG |
---|
673 | if (debugger && debugger->onOptimalPath(*solver)) { |
---|
674 | #if CGL_DEBUG>1 |
---|
675 | const double * optimal = debugger->optimalSolution(); |
---|
676 | CoinPackedVector rpv = thisCut->row(); |
---|
677 | const int n = rpv.getNumElements(); |
---|
678 | const int * indices = rpv.getIndices(); |
---|
679 | const double * elements = rpv.getElements(); |
---|
680 | |
---|
681 | double lb=thisCut->lb(); |
---|
682 | double ub=thisCut->ub(); |
---|
683 | double sum=0.0; |
---|
684 | |
---|
685 | for (int k=0; k<n; k++){ |
---|
686 | int column=indices[k]; |
---|
687 | sum += optimal[column]*elements[k]; |
---|
688 | } |
---|
689 | // is it nearly violated |
---|
690 | if (sum >ub - 1.0e-8 ||sum < lb + 1.0e-8) { |
---|
691 | double violation=CoinMax(sum-ub,lb-sum); |
---|
692 | std::cout<<generatorName_<<" cut with "<<n |
---|
693 | <<" coefficients, nearly cuts off known solutions by "<<violation |
---|
694 | <<", lo="<<lb<<", ub="<<ub<<std::endl; |
---|
695 | for (int k=0; k<n; k++){ |
---|
696 | int column=indices[k]; |
---|
697 | std::cout<<"( "<<column<<" , "<<elements[k]<<" ) "; |
---|
698 | if ((k%4)==3) |
---|
699 | std::cout <<std::endl; |
---|
700 | } |
---|
701 | std::cout <<std::endl; |
---|
702 | std::cout <<"Non zero solution values are"<<std::endl; |
---|
703 | int j=0; |
---|
704 | for (int k=0; k<n; k++){ |
---|
705 | int column=indices[k]; |
---|
706 | if (fabs(optimal[column])>1.0e-9) { |
---|
707 | std::cout<<"( "<<column<<" , "<<optimal[column]<<" ) "; |
---|
708 | if ((j%4)==3) |
---|
709 | std::cout <<std::endl; |
---|
710 | j++; |
---|
711 | } |
---|
712 | } |
---|
713 | std::cout <<std::endl; |
---|
714 | } |
---|
715 | #endif |
---|
716 | assert(!debugger->invalidCut(*thisCut)); |
---|
717 | if(debugger->invalidCut(*thisCut)) |
---|
718 | abort(); |
---|
719 | } |
---|
720 | #endif |
---|
721 | thisCut->mutableRow().setTestForDuplicateIndex(false); |
---|
722 | } |
---|
723 | } |
---|
724 | // Add in saved cuts if violated |
---|
725 | if (false && !depth) { |
---|
726 | const double * solution = solver->getColSolution(); |
---|
727 | double primalTolerance = 1.0e-7; |
---|
728 | int numberCuts = savedCuts_.sizeRowCuts() ; |
---|
729 | for (int k = numberCuts - 1; k >= 0; k--) { |
---|
730 | const OsiRowCut * thisCut = savedCuts_.rowCutPtr(k) ; |
---|
731 | double sum = 0.0; |
---|
732 | int n = thisCut->row().getNumElements(); |
---|
733 | const int * column = thisCut->row().getIndices(); |
---|
734 | const double * element = thisCut->row().getElements(); |
---|
735 | assert (n); |
---|
736 | for (int i = 0; i < n; i++) { |
---|
737 | double value = element[i]; |
---|
738 | sum += value * solution[column[i]]; |
---|
739 | } |
---|
740 | if (sum > thisCut->ub() + primalTolerance) { |
---|
741 | sum = sum - thisCut->ub(); |
---|
742 | } else if (sum < thisCut->lb() - primalTolerance) { |
---|
743 | sum = thisCut->lb() - sum; |
---|
744 | } else { |
---|
745 | sum = 0.0; |
---|
746 | } |
---|
747 | if (sum) { |
---|
748 | // add to candidates and take out here |
---|
749 | cs.insert(*thisCut); |
---|
750 | savedCuts_.eraseRowCut(k); |
---|
751 | } |
---|
752 | } |
---|
753 | } |
---|
754 | if (!atSolution()) { |
---|
755 | int numberRowCutsAfter = cs.sizeRowCuts() ; |
---|
756 | int k ; |
---|
757 | int nEls = 0; |
---|
758 | int nCuts = numberRowCutsAfter - numberRowCutsBefore; |
---|
759 | // Remove NULL cuts! |
---|
760 | int nNull = 0; |
---|
761 | const double * solution = solver->getColSolution(); |
---|
762 | bool feasible = true; |
---|
763 | double primalTolerance = 1.0e-7; |
---|
764 | int shortCut = (depth) ? -1 : generator_->maximumLengthOfCutInTree(); |
---|
765 | for (k = numberRowCutsAfter - 1; k >= numberRowCutsBefore; k--) { |
---|
766 | const OsiRowCut * thisCut = cs.rowCutPtr(k) ; |
---|
767 | double sum = 0.0; |
---|
768 | if (thisCut->lb() <= thisCut->ub()) { |
---|
769 | int n = thisCut->row().getNumElements(); |
---|
770 | if (n <= shortCut) |
---|
771 | numberShortCutsAtRoot_++; |
---|
772 | const int * column = thisCut->row().getIndices(); |
---|
773 | const double * element = thisCut->row().getElements(); |
---|
774 | if (n <= 0) { |
---|
775 | // infeasible cut - give up |
---|
776 | feasible = false; |
---|
777 | break; |
---|
778 | } |
---|
779 | nEls += n; |
---|
780 | for (int i = 0; i < n; i++) { |
---|
781 | double value = element[i]; |
---|
782 | sum += value * solution[column[i]]; |
---|
783 | } |
---|
784 | if (sum > thisCut->ub() + primalTolerance) { |
---|
785 | sum = sum - thisCut->ub(); |
---|
786 | } else if (sum < thisCut->lb() - primalTolerance) { |
---|
787 | sum = thisCut->lb() - sum; |
---|
788 | } else { |
---|
789 | sum = 0.0; |
---|
790 | cs.eraseRowCut(k); |
---|
791 | nNull++; |
---|
792 | } |
---|
793 | } |
---|
794 | } |
---|
795 | //if (nNull) |
---|
796 | //printf("%s has %d cuts and %d elements - %d null!\n",generatorName_, |
---|
797 | // nCuts,nEls,nNull); |
---|
798 | numberRowCutsAfter = cs.sizeRowCuts() ; |
---|
799 | nCuts = numberRowCutsAfter - numberRowCutsBefore; |
---|
800 | nEls = 0; |
---|
801 | for (k = numberRowCutsBefore; k < numberRowCutsAfter; k++) { |
---|
802 | const OsiRowCut * thisCut = cs.rowCutPtr(k) ; |
---|
803 | int n = thisCut->row().getNumElements(); |
---|
804 | nEls += n; |
---|
805 | } |
---|
806 | //printf("%s has %d cuts and %d elements\n",generatorName_, |
---|
807 | // nCuts,nEls); |
---|
808 | int nElsNow = solver->getMatrixByCol()->getNumElements(); |
---|
809 | int numberColumns = solver->getNumCols(); |
---|
810 | int numberRows = solver->getNumRows(); |
---|
811 | //double averagePerRow = static_cast<double>(nElsNow)/ |
---|
812 | //static_cast<double>(numberRows); |
---|
813 | int nAdd; |
---|
814 | int nAdd2; |
---|
815 | int nReasonable; |
---|
816 | if (!model_->parentModel() && depth < 2) { |
---|
817 | if (inaccuracy_ < 3) { |
---|
818 | nAdd = 10000; |
---|
819 | if (pass > 0 && numberColumns > -500) |
---|
820 | nAdd = CoinMin(nAdd, nElsNow + 2 * numberRows); |
---|
821 | } else { |
---|
822 | nAdd = 10000; |
---|
823 | if (pass > 0) |
---|
824 | nAdd = CoinMin(nAdd, nElsNow + 2 * numberRows); |
---|
825 | } |
---|
826 | nAdd2 = 5 * numberColumns; |
---|
827 | nReasonable = CoinMax(nAdd2, nElsNow / 8 + nAdd); |
---|
828 | if (!depth && !pass) { |
---|
829 | // allow more |
---|
830 | nAdd += nElsNow / 2; |
---|
831 | nAdd2 += nElsNow / 2; |
---|
832 | nReasonable += nElsNow / 2; |
---|
833 | } |
---|
834 | //if (!depth&&ineffectualCuts()) |
---|
835 | //nReasonable *= 2; |
---|
836 | } else { |
---|
837 | nAdd = 200; |
---|
838 | nAdd2 = 2 * numberColumns; |
---|
839 | nReasonable = CoinMax(nAdd2, nElsNow / 8 + nAdd); |
---|
840 | } |
---|
841 | //#define UNS_WEIGHT 0.1 |
---|
842 | #ifdef UNS_WEIGHT |
---|
843 | const double * colLower = solver->getColLower(); |
---|
844 | const double * colUpper = solver->getColUpper(); |
---|
845 | #endif |
---|
846 | if (/*nEls>CoinMax(nAdd2,nElsNow/8+nAdd)*/nCuts && feasible) { |
---|
847 | //printf("need to remove cuts\n"); |
---|
848 | // just add most effective |
---|
849 | #ifndef JJF_ONE |
---|
850 | int nDelete = nEls - nReasonable; |
---|
851 | |
---|
852 | nElsNow = nEls; |
---|
853 | double * sort = new double [nCuts]; |
---|
854 | int * which = new int [nCuts]; |
---|
855 | // For parallel cuts |
---|
856 | double * element2 = new double [numberColumns]; |
---|
857 | //#define USE_OBJECTIVE 2 |
---|
858 | #ifdef USE_OBJECTIVE |
---|
859 | const double *objective = solver->getObjCoefficients() ; |
---|
860 | #if USE_OBJECTIVE>1 |
---|
861 | double objNorm = 0.0; |
---|
862 | for (int i = 0; i < numberColumns; i++) |
---|
863 | objNorm += objective[i] * objective[i]; |
---|
864 | if (objNorm) |
---|
865 | objNorm = 1.0 / sqrt(objNorm); |
---|
866 | else |
---|
867 | objNorm = 1.0; |
---|
868 | objNorm *= 0.01; // downgrade |
---|
869 | #endif |
---|
870 | #endif |
---|
871 | CoinZeroN(element2, numberColumns); |
---|
872 | for (k = numberRowCutsBefore; k < numberRowCutsAfter; k++) { |
---|
873 | const OsiRowCut * thisCut = cs.rowCutPtr(k) ; |
---|
874 | double sum = 0.0; |
---|
875 | if (thisCut->lb() <= thisCut->ub()) { |
---|
876 | int n = thisCut->row().getNumElements(); |
---|
877 | const int * column = thisCut->row().getIndices(); |
---|
878 | const double * element = thisCut->row().getElements(); |
---|
879 | assert (n); |
---|
880 | #ifdef UNS_WEIGHT |
---|
881 | double normU = 0.0; |
---|
882 | double norm = 1.0e-3; |
---|
883 | int nU = 0; |
---|
884 | for (int i = 0; i < n; i++) { |
---|
885 | double value = element[i]; |
---|
886 | int iColumn = column[i]; |
---|
887 | double solValue = solution[iColumn]; |
---|
888 | sum += value * solValue; |
---|
889 | value *= value; |
---|
890 | norm += value; |
---|
891 | if (solValue > colLower[iColumn] + 1.0e-6 && |
---|
892 | solValue < colUpper[iColumn] - 1.0e-6) { |
---|
893 | normU += value; |
---|
894 | nU++; |
---|
895 | } |
---|
896 | } |
---|
897 | #ifdef JJF_ZERO |
---|
898 | int nS = n - nU; |
---|
899 | if (numberColumns > 20000) { |
---|
900 | if (nS > 50) { |
---|
901 | double ratio = 50.0 / nS; |
---|
902 | normU /= ratio; |
---|
903 | } |
---|
904 | } |
---|
905 | #endif |
---|
906 | norm += UNS_WEIGHT * (normU - norm); |
---|
907 | #else |
---|
908 | double norm = 1.0e-3; |
---|
909 | #ifdef USE_OBJECTIVE |
---|
910 | double obj = 0.0; |
---|
911 | #endif |
---|
912 | for (int i = 0; i < n; i++) { |
---|
913 | int iColumn = column[i]; |
---|
914 | double value = element[i]; |
---|
915 | sum += value * solution[iColumn]; |
---|
916 | norm += value * value; |
---|
917 | #ifdef USE_OBJECTIVE |
---|
918 | obj += value * objective[iColumn]; |
---|
919 | #endif |
---|
920 | } |
---|
921 | #endif |
---|
922 | if (sum > thisCut->ub()) { |
---|
923 | sum = sum - thisCut->ub(); |
---|
924 | } else if (sum < thisCut->lb()) { |
---|
925 | sum = thisCut->lb() - sum; |
---|
926 | } else { |
---|
927 | sum = 0.0; |
---|
928 | } |
---|
929 | #ifdef USE_OBJECTIVE |
---|
930 | if (sum) { |
---|
931 | #if USE_OBJECTIVE==1 |
---|
932 | obj = CoinMax(1.0e-6, fabs(obj)); |
---|
933 | norm = sqrt(obj * norm); |
---|
934 | //sum += fabs(obj)*invObjNorm; |
---|
935 | //printf("sum %g norm %g normobj %g invNorm %g mod %g\n", |
---|
936 | // sum,norm,obj,invObjNorm,obj*invObjNorm); |
---|
937 | // normalize |
---|
938 | sum /= sqrt(norm); |
---|
939 | #else |
---|
940 | // normalize |
---|
941 | norm = 1.0 / sqrt(norm); |
---|
942 | sum = (sum + objNorm * obj) * norm; |
---|
943 | #endif |
---|
944 | } |
---|
945 | #else |
---|
946 | // normalize |
---|
947 | sum /= sqrt(norm); |
---|
948 | #endif |
---|
949 | //sum /= pow(norm,0.3); |
---|
950 | // adjust for length |
---|
951 | //sum /= pow(reinterpret_cast<double>(n),0.2); |
---|
952 | //sum /= sqrt((double) n); |
---|
953 | // randomize |
---|
954 | //double randomNumber = |
---|
955 | //model_->randomNumberGenerator()->randomDouble(); |
---|
956 | //sum *= (0.5+randomNumber); |
---|
957 | } else { |
---|
958 | // keep |
---|
959 | sum = COIN_DBL_MAX; |
---|
960 | } |
---|
961 | sort[k-numberRowCutsBefore] = sum; |
---|
962 | which[k-numberRowCutsBefore] = k; |
---|
963 | } |
---|
964 | CoinSort_2(sort, sort + nCuts, which); |
---|
965 | // Now see which ones are too similar |
---|
966 | int nParallel = 0; |
---|
967 | double testValue = (depth > 1) ? 0.99 : 0.999999; |
---|
968 | for (k = 0; k < nCuts; k++) { |
---|
969 | int j = which[k]; |
---|
970 | const OsiRowCut * thisCut = cs.rowCutPtr(j) ; |
---|
971 | if (thisCut->lb() > thisCut->ub()) |
---|
972 | break; // cut is infeasible |
---|
973 | int n = thisCut->row().getNumElements(); |
---|
974 | const int * column = thisCut->row().getIndices(); |
---|
975 | const double * element = thisCut->row().getElements(); |
---|
976 | assert (n); |
---|
977 | double norm = 0.0; |
---|
978 | double lb = thisCut->lb(); |
---|
979 | double ub = thisCut->ub(); |
---|
980 | for (int i = 0; i < n; i++) { |
---|
981 | double value = element[i]; |
---|
982 | element2[column[i]] = value; |
---|
983 | norm += value * value; |
---|
984 | } |
---|
985 | int kkk = CoinMin(nCuts, k + 5); |
---|
986 | for (int kk = k + 1; kk < kkk; kk++) { |
---|
987 | int jj = which[kk]; |
---|
988 | const OsiRowCut * thisCut2 = cs.rowCutPtr(jj) ; |
---|
989 | if (thisCut2->lb() > thisCut2->ub()) |
---|
990 | break; // cut is infeasible |
---|
991 | int nB = thisCut2->row().getNumElements(); |
---|
992 | const int * columnB = thisCut2->row().getIndices(); |
---|
993 | const double * elementB = thisCut2->row().getElements(); |
---|
994 | assert (nB); |
---|
995 | double normB = 0.0; |
---|
996 | double product = 0.0; |
---|
997 | for (int i = 0; i < nB; i++) { |
---|
998 | double value = elementB[i]; |
---|
999 | normB += value * value; |
---|
1000 | product += value * element2[columnB[i]]; |
---|
1001 | } |
---|
1002 | if (product > 0.0 && product*product > testValue*norm*normB) { |
---|
1003 | bool parallel = true; |
---|
1004 | double lbB = thisCut2->lb(); |
---|
1005 | double ubB = thisCut2->ub(); |
---|
1006 | if ((lb < -1.0e20 && lbB > -1.0e20) || |
---|
1007 | (lbB < -1.0e20 && lb > -1.0e20)) |
---|
1008 | parallel = false; |
---|
1009 | double tolerance; |
---|
1010 | tolerance = CoinMax(fabs(lb), fabs(lbB)) + 1.0e-6; |
---|
1011 | if (fabs(lb - lbB) > tolerance) |
---|
1012 | parallel = false; |
---|
1013 | if ((ub > 1.0e20 && ubB < 1.0e20) || |
---|
1014 | (ubB > 1.0e20 && ub < 1.0e20)) |
---|
1015 | parallel = false; |
---|
1016 | tolerance = CoinMax(fabs(ub), fabs(ubB)) + 1.0e-6; |
---|
1017 | if (fabs(ub - ubB) > tolerance) |
---|
1018 | parallel = false; |
---|
1019 | if (parallel) { |
---|
1020 | nParallel++; |
---|
1021 | sort[k] = 0.0; |
---|
1022 | break; |
---|
1023 | } |
---|
1024 | } |
---|
1025 | } |
---|
1026 | for (int i = 0; i < n; i++) { |
---|
1027 | element2[column[i]] = 0.0; |
---|
1028 | } |
---|
1029 | } |
---|
1030 | delete [] element2; |
---|
1031 | CoinSort_2(sort, sort + nCuts, which); |
---|
1032 | k = 0; |
---|
1033 | while (nDelete > 0 || !sort[k]) { |
---|
1034 | int iCut = which[k]; |
---|
1035 | const OsiRowCut * thisCut = cs.rowCutPtr(iCut) ; |
---|
1036 | int n = thisCut->row().getNumElements(); |
---|
1037 | // may be best, just to save if short |
---|
1038 | if (false && n && sort[k]) { |
---|
1039 | // add to saved cuts |
---|
1040 | savedCuts_.insert(*thisCut); |
---|
1041 | } |
---|
1042 | nDelete -= n; |
---|
1043 | k++; |
---|
1044 | if (k >= nCuts) |
---|
1045 | break; |
---|
1046 | } |
---|
1047 | std::sort(which, which + k); |
---|
1048 | k--; |
---|
1049 | for (; k >= 0; k--) { |
---|
1050 | cs.eraseRowCut(which[k]); |
---|
1051 | } |
---|
1052 | delete [] sort; |
---|
1053 | delete [] which; |
---|
1054 | numberRowCutsAfter = cs.sizeRowCuts() ; |
---|
1055 | #else |
---|
1056 | double * norm = new double [nCuts]; |
---|
1057 | int * which = new int [2*nCuts]; |
---|
1058 | double * score = new double [nCuts]; |
---|
1059 | double * ortho = new double [nCuts]; |
---|
1060 | int nIn = 0; |
---|
1061 | int nOut = nCuts; |
---|
1062 | // For parallel cuts |
---|
1063 | double * element2 = new double [numberColumns]; |
---|
1064 | const double *objective = solver->getObjCoefficients() ; |
---|
1065 | double objNorm = 0.0; |
---|
1066 | for (int i = 0; i < numberColumns; i++) |
---|
1067 | objNorm += objective[i] * objective[i]; |
---|
1068 | if (objNorm) |
---|
1069 | objNorm = 1.0 / sqrt(objNorm); |
---|
1070 | else |
---|
1071 | objNorm = 1.0; |
---|
1072 | objNorm *= 0.1; // weight of 0.1 |
---|
1073 | CoinZeroN(element2, numberColumns); |
---|
1074 | int numberRowCuts = numberRowCutsAfter - numberRowCutsBefore; |
---|
1075 | int iBest = -1; |
---|
1076 | double best = 0.0; |
---|
1077 | int nPossible = 0; |
---|
1078 | double testValue = (depth > 1) ? 0.7 : 0.5; |
---|
1079 | for (k = 0; k < numberRowCuts; k++) { |
---|
1080 | const OsiRowCut * thisCut = cs.rowCutPtr(k + numberRowCutsBefore) ; |
---|
1081 | double sum = 0.0; |
---|
1082 | if (thisCut->lb() <= thisCut->ub()) { |
---|
1083 | int n = thisCut->row().getNumElements(); |
---|
1084 | const int * column = thisCut->row().getIndices(); |
---|
1085 | const double * element = thisCut->row().getElements(); |
---|
1086 | assert (n); |
---|
1087 | double normThis = 1.0e-6; |
---|
1088 | double obj = 0.0; |
---|
1089 | for (int i = 0; i < n; i++) { |
---|
1090 | int iColumn = column[i]; |
---|
1091 | double value = element[i]; |
---|
1092 | sum += value * solution[iColumn]; |
---|
1093 | normThis += value * value; |
---|
1094 | obj += value * objective[iColumn]; |
---|
1095 | } |
---|
1096 | if (sum > thisCut->ub()) { |
---|
1097 | sum = sum - thisCut->ub(); |
---|
1098 | } else if (sum < thisCut->lb()) { |
---|
1099 | sum = thisCut->lb() - sum; |
---|
1100 | } else { |
---|
1101 | sum = 0.0; |
---|
1102 | } |
---|
1103 | if (sum) { |
---|
1104 | normThis = 1.0 / sqrt(normThis); |
---|
1105 | norm[k] = normThis; |
---|
1106 | sum *= normThis; |
---|
1107 | obj *= normThis; |
---|
1108 | score[k] = sum + obj * objNorm; |
---|
1109 | ortho[k] = 1.0; |
---|
1110 | } |
---|
1111 | } else { |
---|
1112 | // keep and discard others |
---|
1113 | nIn = 1; |
---|
1114 | which[0] = k; |
---|
1115 | for (int j = 0; j < numberRowCuts; j++) { |
---|
1116 | if (j != k) |
---|
1117 | which[nOut++] = j; |
---|
1118 | } |
---|
1119 | iBest = -1; |
---|
1120 | break; |
---|
1121 | } |
---|
1122 | if (sum) { |
---|
1123 | if (score[k] > best) { |
---|
1124 | best = score[k]; |
---|
1125 | iBest = nPossible; |
---|
1126 | } |
---|
1127 | which[nPossible++] = k; |
---|
1128 | } else { |
---|
1129 | which[nOut++] = k; |
---|
1130 | } |
---|
1131 | } |
---|
1132 | while (iBest >= 0) { |
---|
1133 | int kBest = which[iBest]; |
---|
1134 | int j = which[nIn]; |
---|
1135 | which[iBest] = j; |
---|
1136 | which[nIn++] = kBest; |
---|
1137 | const OsiRowCut * thisCut = cs.rowCutPtr(kBest + numberRowCutsBefore) ; |
---|
1138 | int n = thisCut->row().getNumElements(); |
---|
1139 | nReasonable -= n; |
---|
1140 | if (nReasonable <= 0) { |
---|
1141 | for (k = nIn; k < nPossible; k++) |
---|
1142 | which[nOut++] = which[k]; |
---|
1143 | break; |
---|
1144 | } |
---|
1145 | // Now see which ones are too similar and choose next |
---|
1146 | iBest = -1; |
---|
1147 | best = 0.0; |
---|
1148 | int nOld = nPossible; |
---|
1149 | nPossible = nIn; |
---|
1150 | const int * column = thisCut->row().getIndices(); |
---|
1151 | const double * element = thisCut->row().getElements(); |
---|
1152 | assert (n); |
---|
1153 | double normNew = norm[kBest]; |
---|
1154 | for (int i = 0; i < n; i++) { |
---|
1155 | double value = element[i]; |
---|
1156 | element2[column[i]] = value; |
---|
1157 | } |
---|
1158 | for (int j = nIn; j < nOld; j++) { |
---|
1159 | k = which[j]; |
---|
1160 | const OsiRowCut * thisCut2 = cs.rowCutPtr(k + numberRowCutsBefore) ; |
---|
1161 | int nB = thisCut2->row().getNumElements(); |
---|
1162 | const int * columnB = thisCut2->row().getIndices(); |
---|
1163 | const double * elementB = thisCut2->row().getElements(); |
---|
1164 | assert (nB); |
---|
1165 | double normB = norm[k]; |
---|
1166 | double product = 0.0; |
---|
1167 | for (int i = 0; i < nB; i++) { |
---|
1168 | double value = elementB[i]; |
---|
1169 | product += value * element2[columnB[i]]; |
---|
1170 | } |
---|
1171 | double orthoScore = 1.0 - product * normNew * normB; |
---|
1172 | if (orthoScore >= testValue) { |
---|
1173 | ortho[k] = CoinMin(orthoScore, ortho[k]); |
---|
1174 | double test = score[k] + ortho[k]; |
---|
1175 | if (test > best) { |
---|
1176 | best = score[k]; |
---|
1177 | iBest = nPossible; |
---|
1178 | } |
---|
1179 | which[nPossible++] = k; |
---|
1180 | } else { |
---|
1181 | which[nOut++] = k; |
---|
1182 | } |
---|
1183 | } |
---|
1184 | for (int i = 0; i < n; i++) { |
---|
1185 | element2[column[i]] = 0.0; |
---|
1186 | } |
---|
1187 | } |
---|
1188 | delete [] score; |
---|
1189 | delete [] ortho; |
---|
1190 | std::sort(which + nCuts, which + nOut); |
---|
1191 | k = nOut - 1; |
---|
1192 | for (; k >= nCuts; k--) { |
---|
1193 | cs.eraseRowCut(which[k] + numberRowCutsBefore); |
---|
1194 | } |
---|
1195 | delete [] norm; |
---|
1196 | delete [] which; |
---|
1197 | numberRowCutsAfter = cs.sizeRowCuts() ; |
---|
1198 | #endif |
---|
1199 | } |
---|
1200 | } |
---|
1201 | #ifdef CBC_DEBUG |
---|
1202 | { |
---|
1203 | int numberRowCutsAfter = cs.sizeRowCuts() ; |
---|
1204 | int k ; |
---|
1205 | int nBad = 0; |
---|
1206 | for (k = numberRowCutsBefore; k < numberRowCutsAfter; k++) { |
---|
1207 | OsiRowCut thisCut = cs.rowCut(k) ; |
---|
1208 | if (thisCut.lb() > thisCut.ub() || |
---|
1209 | thisCut.lb() > 1.0e8 || |
---|
1210 | thisCut.ub() < -1.0e8) |
---|
1211 | printf("cut from %s has bounds %g and %g!\n", |
---|
1212 | generatorName_, thisCut.lb(), thisCut.ub()); |
---|
1213 | if (thisCut.lb() <= thisCut.ub()) { |
---|
1214 | /* check size of elements. |
---|
1215 | We can allow smaller but this helps debug generators as it |
---|
1216 | is unsafe to have small elements */ |
---|
1217 | int n = thisCut.row().getNumElements(); |
---|
1218 | const int * column = thisCut.row().getIndices(); |
---|
1219 | const double * element = thisCut.row().getElements(); |
---|
1220 | assert (n); |
---|
1221 | for (int i = 0; i < n; i++) { |
---|
1222 | double value = element[i]; |
---|
1223 | if (fabs(value) <= 1.0e-12 || fabs(value) >= 1.0e20) |
---|
1224 | nBad++; |
---|
1225 | } |
---|
1226 | } |
---|
1227 | if (nBad) |
---|
1228 | printf("Cut generator %s produced %d cuts of which %d had tiny or large elements\n", |
---|
1229 | generatorName_, numberRowCutsAfter - numberRowCutsBefore, nBad); |
---|
1230 | } |
---|
1231 | } |
---|
1232 | #endif |
---|
1233 | int numberRowCutsAfter = cs.sizeRowCuts() ; |
---|
1234 | int numberColumnCutsAfter = cs.sizeColCuts() ; |
---|
1235 | if (numberRowCutsBefore < numberRowCutsAfter) { |
---|
1236 | for (int k = numberRowCutsBefore; k < numberRowCutsAfter; k++) { |
---|
1237 | OsiRowCut thisCut = cs.rowCut(k) ; |
---|
1238 | int n = thisCut.row().getNumElements(); |
---|
1239 | numberElements_ += n; |
---|
1240 | } |
---|
1241 | #ifdef JJF_ZERO |
---|
1242 | printf("generator %s generated %d row cuts\n", |
---|
1243 | generatorName_, numberRowCutsAfter - numberRowCutsBefore); |
---|
1244 | #endif |
---|
1245 | numberCuts_ += numberRowCutsAfter - numberRowCutsBefore; |
---|
1246 | } |
---|
1247 | if (numberColumnCutsBefore < numberColumnCutsAfter) { |
---|
1248 | #ifdef JJF_ZERO |
---|
1249 | printf("generator %s generated %d column cuts\n", |
---|
1250 | generatorName_, numberColumnCutsAfter - numberColumnCutsBefore); |
---|
1251 | #endif |
---|
1252 | numberColumnCuts_ += numberColumnCutsAfter - numberColumnCutsBefore; |
---|
1253 | } |
---|
1254 | if (timing()) |
---|
1255 | timeInCutGenerator_ += CoinCpuTime() - time1; |
---|
1256 | // switch off if first time and no good |
---|
1257 | if (node == NULL && !pass ) { |
---|
1258 | if (numberRowCutsAfter - numberRowCutsBefore |
---|
1259 | < switchOffIfLessThan_ /*&& numberCuts_ < switchOffIfLessThan_*/) { |
---|
1260 | // switch off |
---|
1261 | maximumTries_ = 0; |
---|
1262 | whenCutGenerator_=-100; |
---|
1263 | //whenCutGenerator_ = -100; |
---|
1264 | //whenCutGeneratorInSub_ = -200; |
---|
1265 | } |
---|
1266 | } |
---|
1267 | if (maximumTries_>0) { |
---|
1268 | maximumTries_--; |
---|
1269 | if (!maximumTries_) |
---|
1270 | whenCutGenerator_=-100; |
---|
1271 | } |
---|
1272 | } |
---|
1273 | return returnCode; |
---|
1274 | } |
---|
1275 | void |
---|
1276 | CbcCutGenerator::setHowOften(int howOften) |
---|
1277 | { |
---|
1278 | |
---|
1279 | if (howOften >= 1000000) { |
---|
1280 | // leave Probing every SCANCUTS_PROBING |
---|
1281 | howOften = howOften % 1000000; |
---|
1282 | CglProbing* generator = |
---|
1283 | dynamic_cast<CglProbing*>(generator_); |
---|
1284 | |
---|
1285 | if (generator && howOften > SCANCUTS_PROBING) |
---|
1286 | howOften = SCANCUTS_PROBING + 1000000; |
---|
1287 | else |
---|
1288 | howOften += 1000000; |
---|
1289 | } |
---|
1290 | whenCutGenerator_ = howOften; |
---|
1291 | } |
---|
1292 | void |
---|
1293 | CbcCutGenerator::setWhatDepth(int value) |
---|
1294 | { |
---|
1295 | depthCutGenerator_ = value; |
---|
1296 | } |
---|
1297 | void |
---|
1298 | CbcCutGenerator::setWhatDepthInSub(int value) |
---|
1299 | { |
---|
1300 | depthCutGeneratorInSub_ = value; |
---|
1301 | } |
---|
1302 | // Add in statistics from other |
---|
1303 | void |
---|
1304 | CbcCutGenerator::addStatistics(const CbcCutGenerator * other) |
---|
1305 | { |
---|
1306 | // Time in cut generator |
---|
1307 | timeInCutGenerator_ += other->timeInCutGenerator_; |
---|
1308 | // Number times cut generator entered |
---|
1309 | numberTimes_ += other->numberTimes_; |
---|
1310 | // Total number of cuts added |
---|
1311 | numberCuts_ += other->numberCuts_; |
---|
1312 | // Total number of elements added |
---|
1313 | numberElements_ += other->numberElements_; |
---|
1314 | // Total number of column cuts added |
---|
1315 | numberColumnCuts_ += other->numberColumnCuts_; |
---|
1316 | // Total number of cuts active after (at end of n cut passes at each node) |
---|
1317 | numberCutsActive_ += other->numberCutsActive_; |
---|
1318 | // Number of cuts generated at root |
---|
1319 | numberCutsAtRoot_ += other->numberCutsAtRoot_; |
---|
1320 | // Number of cuts active at root |
---|
1321 | numberActiveCutsAtRoot_ += other->numberActiveCutsAtRoot_; |
---|
1322 | // Number of short cuts at root |
---|
1323 | numberShortCutsAtRoot_ += other->numberShortCutsAtRoot_; |
---|
1324 | } |
---|
1325 | // Scale back statistics by factor |
---|
1326 | void |
---|
1327 | CbcCutGenerator::scaleBackStatistics(int factor) |
---|
1328 | { |
---|
1329 | // leave time |
---|
1330 | // Number times cut generator entered |
---|
1331 | numberTimes_ = (numberTimes_+factor-1)/factor; |
---|
1332 | // Total number of cuts added |
---|
1333 | numberCuts_ = (numberCuts_+factor-1)/factor; |
---|
1334 | // Total number of elements added |
---|
1335 | numberElements_ = (numberElements_+factor-1)/factor; |
---|
1336 | // Total number of column cuts added |
---|
1337 | numberColumnCuts_ = (numberColumnCuts_+factor-1)/factor; |
---|
1338 | // Total number of cuts active after (at end of n cut passes at each node) |
---|
1339 | numberCutsActive_ = (numberCutsActive_+factor-1)/factor; |
---|
1340 | // Number of cuts generated at root |
---|
1341 | numberCutsAtRoot_ = (numberCutsAtRoot_+factor-1)/factor; |
---|
1342 | // Number of cuts active at root |
---|
1343 | numberActiveCutsAtRoot_ = (numberActiveCutsAtRoot_+factor-1)/factor; |
---|
1344 | // Number of short cuts at root |
---|
1345 | numberShortCutsAtRoot_ = (numberShortCutsAtRoot_+factor-1)/factor; |
---|
1346 | } |
---|
1347 | // Create C++ lines to get to current state |
---|
1348 | void |
---|
1349 | CbcCutGenerator::generateTuning( FILE * fp) |
---|
1350 | { |
---|
1351 | fprintf(fp, "// Cbc tuning for generator %s\n", generatorName_); |
---|
1352 | fprintf(fp, " generator->setHowOften(%d);\n", whenCutGenerator_); |
---|
1353 | fprintf(fp, " generator->setSwitchOffIfLessThan(%d);\n", switchOffIfLessThan_); |
---|
1354 | fprintf(fp, " generator->setWhatDepth(%d);\n", depthCutGenerator_); |
---|
1355 | fprintf(fp, " generator->setInaccuracy(%d);\n", inaccuracy_); |
---|
1356 | if (timing()) |
---|
1357 | fprintf(fp, " generator->setTiming(true);\n"); |
---|
1358 | if (normal()) |
---|
1359 | fprintf(fp, " generator->setNormal(true);\n"); |
---|
1360 | if (atSolution()) |
---|
1361 | fprintf(fp, " generator->setAtSolution(true);\n"); |
---|
1362 | if (whenInfeasible()) |
---|
1363 | fprintf(fp, " generator->setWhenInfeasible(true);\n"); |
---|
1364 | if (needsOptimalBasis()) |
---|
1365 | fprintf(fp, " generator->setNeedsOptimalBasis(true);\n"); |
---|
1366 | if (mustCallAgain()) |
---|
1367 | fprintf(fp, " generator->setMustCallAgain(true);\n"); |
---|
1368 | if (whetherToUse()) |
---|
1369 | fprintf(fp, " generator->setWhetherToUse(true);\n"); |
---|
1370 | } |
---|
1371 | |
---|
1372 | |
---|
1373 | |
---|