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