1 | /* $Id: CbcHeuristic.cpp 2280 2016-06-14 14:39:54Z forrest $ */ |
---|
2 | // Copyright (C) 2002, 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 | |
---|
11 | #include "CbcConfig.h" |
---|
12 | |
---|
13 | #include <cassert> |
---|
14 | #include <cstdlib> |
---|
15 | #include <cmath> |
---|
16 | #include <cfloat> |
---|
17 | |
---|
18 | //#define PRINT_DEBUG |
---|
19 | #ifdef COIN_HAS_CLP |
---|
20 | #include "OsiClpSolverInterface.hpp" |
---|
21 | #endif |
---|
22 | #include "CbcModel.hpp" |
---|
23 | #include "CbcMessage.hpp" |
---|
24 | #include "CbcHeuristic.hpp" |
---|
25 | #include "CbcHeuristicFPump.hpp" |
---|
26 | #include "CbcHeuristicRINS.hpp" |
---|
27 | #include "CbcEventHandler.hpp" |
---|
28 | #include "CbcStrategy.hpp" |
---|
29 | #include "CglPreProcess.hpp" |
---|
30 | #include "CglGomory.hpp" |
---|
31 | #include "CglProbing.hpp" |
---|
32 | #include "OsiAuxInfo.hpp" |
---|
33 | #include "OsiRowCutDebugger.hpp" |
---|
34 | #include "OsiPresolve.hpp" |
---|
35 | #include "CbcBranchActual.hpp" |
---|
36 | #include "CbcCutGenerator.hpp" |
---|
37 | //============================================================================== |
---|
38 | |
---|
39 | CbcHeuristicNode::CbcHeuristicNode(const CbcHeuristicNode& rhs) |
---|
40 | { |
---|
41 | numObjects_ = rhs.numObjects_; |
---|
42 | brObj_ = new CbcBranchingObject*[numObjects_]; |
---|
43 | for (int i = 0; i < numObjects_; ++i) { |
---|
44 | brObj_[i] = rhs.brObj_[i]->clone(); |
---|
45 | } |
---|
46 | } |
---|
47 | |
---|
48 | void |
---|
49 | CbcHeuristicNodeList::gutsOfDelete() |
---|
50 | { |
---|
51 | for (int i = (static_cast<int>(nodes_.size())) - 1; i >= 0; --i) { |
---|
52 | delete nodes_[i]; |
---|
53 | } |
---|
54 | } |
---|
55 | |
---|
56 | void |
---|
57 | CbcHeuristicNodeList::gutsOfCopy(const CbcHeuristicNodeList& rhs) |
---|
58 | { |
---|
59 | append(rhs); |
---|
60 | } |
---|
61 | |
---|
62 | CbcHeuristicNodeList::CbcHeuristicNodeList(const CbcHeuristicNodeList& rhs) |
---|
63 | { |
---|
64 | gutsOfCopy(rhs); |
---|
65 | } |
---|
66 | |
---|
67 | CbcHeuristicNodeList& CbcHeuristicNodeList::operator= |
---|
68 | (const CbcHeuristicNodeList & rhs) |
---|
69 | { |
---|
70 | if (this != &rhs) { |
---|
71 | gutsOfDelete(); |
---|
72 | gutsOfCopy(rhs); |
---|
73 | } |
---|
74 | return *this; |
---|
75 | } |
---|
76 | |
---|
77 | CbcHeuristicNodeList::~CbcHeuristicNodeList() |
---|
78 | { |
---|
79 | gutsOfDelete(); |
---|
80 | } |
---|
81 | |
---|
82 | void |
---|
83 | CbcHeuristicNodeList::append(CbcHeuristicNode*& node) |
---|
84 | { |
---|
85 | nodes_.push_back(node); |
---|
86 | node = NULL; |
---|
87 | } |
---|
88 | |
---|
89 | void |
---|
90 | CbcHeuristicNodeList::append(const CbcHeuristicNodeList& nodes) |
---|
91 | { |
---|
92 | nodes_.reserve(nodes_.size() + nodes.size()); |
---|
93 | for (int i = 0; i < nodes.size(); ++i) { |
---|
94 | CbcHeuristicNode* node = new CbcHeuristicNode(*nodes.node(i)); |
---|
95 | append(node); |
---|
96 | } |
---|
97 | } |
---|
98 | |
---|
99 | //============================================================================== |
---|
100 | #define DEFAULT_WHERE ((255-2-16)*(1+256)) |
---|
101 | // Default Constructor |
---|
102 | CbcHeuristic::CbcHeuristic() : |
---|
103 | model_(NULL), |
---|
104 | when_(2), |
---|
105 | numberNodes_(200), |
---|
106 | feasibilityPumpOptions_(-1), |
---|
107 | fractionSmall_(1.0), |
---|
108 | heuristicName_("Unknown"), |
---|
109 | howOften_(1), |
---|
110 | decayFactor_(0.0), |
---|
111 | switches_(0), |
---|
112 | whereFrom_(DEFAULT_WHERE), |
---|
113 | shallowDepth_(1), |
---|
114 | howOftenShallow_(1), |
---|
115 | numInvocationsInShallow_(0), |
---|
116 | numInvocationsInDeep_(0), |
---|
117 | lastRunDeep_(0), |
---|
118 | numRuns_(0), |
---|
119 | minDistanceToRun_(1), |
---|
120 | runNodes_(), |
---|
121 | numCouldRun_(0), |
---|
122 | numberSolutionsFound_(0), |
---|
123 | numberNodesDone_(0), |
---|
124 | inputSolution_(NULL) |
---|
125 | { |
---|
126 | // As CbcHeuristic virtual need to modify .cpp if above change |
---|
127 | } |
---|
128 | |
---|
129 | // Constructor from model |
---|
130 | CbcHeuristic::CbcHeuristic(CbcModel & model) : |
---|
131 | model_(&model), |
---|
132 | when_(2), |
---|
133 | numberNodes_(200), |
---|
134 | feasibilityPumpOptions_(-1), |
---|
135 | fractionSmall_(1.0), |
---|
136 | heuristicName_("Unknown"), |
---|
137 | howOften_(1), |
---|
138 | decayFactor_(0.0), |
---|
139 | switches_(0), |
---|
140 | whereFrom_(DEFAULT_WHERE), |
---|
141 | shallowDepth_(1), |
---|
142 | howOftenShallow_(1), |
---|
143 | numInvocationsInShallow_(0), |
---|
144 | numInvocationsInDeep_(0), |
---|
145 | lastRunDeep_(0), |
---|
146 | numRuns_(0), |
---|
147 | minDistanceToRun_(1), |
---|
148 | runNodes_(), |
---|
149 | numCouldRun_(0), |
---|
150 | numberSolutionsFound_(0), |
---|
151 | numberNodesDone_(0), |
---|
152 | inputSolution_(NULL) |
---|
153 | { |
---|
154 | } |
---|
155 | |
---|
156 | void |
---|
157 | CbcHeuristic::gutsOfCopy(const CbcHeuristic & rhs) |
---|
158 | { |
---|
159 | model_ = rhs.model_; |
---|
160 | when_ = rhs.when_; |
---|
161 | numberNodes_ = rhs.numberNodes_; |
---|
162 | feasibilityPumpOptions_ = rhs.feasibilityPumpOptions_; |
---|
163 | fractionSmall_ = rhs.fractionSmall_; |
---|
164 | randomNumberGenerator_ = rhs.randomNumberGenerator_; |
---|
165 | heuristicName_ = rhs.heuristicName_; |
---|
166 | howOften_ = rhs.howOften_; |
---|
167 | decayFactor_ = rhs.decayFactor_; |
---|
168 | switches_ = rhs.switches_; |
---|
169 | whereFrom_ = rhs.whereFrom_; |
---|
170 | shallowDepth_ = rhs.shallowDepth_; |
---|
171 | howOftenShallow_ = rhs.howOftenShallow_; |
---|
172 | numInvocationsInShallow_ = rhs.numInvocationsInShallow_; |
---|
173 | numInvocationsInDeep_ = rhs.numInvocationsInDeep_; |
---|
174 | lastRunDeep_ = rhs.lastRunDeep_; |
---|
175 | numRuns_ = rhs.numRuns_; |
---|
176 | numCouldRun_ = rhs.numCouldRun_; |
---|
177 | minDistanceToRun_ = rhs.minDistanceToRun_; |
---|
178 | runNodes_ = rhs.runNodes_; |
---|
179 | numberSolutionsFound_ = rhs.numberSolutionsFound_; |
---|
180 | numberNodesDone_ = rhs.numberNodesDone_; |
---|
181 | if (rhs.inputSolution_) { |
---|
182 | int numberColumns = model_->getNumCols(); |
---|
183 | setInputSolution(rhs.inputSolution_, rhs.inputSolution_[numberColumns]); |
---|
184 | } |
---|
185 | } |
---|
186 | // Copy constructor |
---|
187 | CbcHeuristic::CbcHeuristic(const CbcHeuristic & rhs) |
---|
188 | { |
---|
189 | inputSolution_ = NULL; |
---|
190 | gutsOfCopy(rhs); |
---|
191 | } |
---|
192 | |
---|
193 | // Assignment operator |
---|
194 | CbcHeuristic & |
---|
195 | CbcHeuristic::operator=( const CbcHeuristic & rhs) |
---|
196 | { |
---|
197 | if (this != &rhs) { |
---|
198 | gutsOfDelete(); |
---|
199 | gutsOfCopy(rhs); |
---|
200 | } |
---|
201 | return *this; |
---|
202 | } |
---|
203 | |
---|
204 | void CbcHeurDebugNodes(CbcModel* model_) |
---|
205 | { |
---|
206 | CbcNode* node = model_->currentNode(); |
---|
207 | CbcNodeInfo* nodeInfo = node->nodeInfo(); |
---|
208 | std::cout << "===============================================================\n"; |
---|
209 | while (nodeInfo) { |
---|
210 | const CbcNode* node = nodeInfo->owner(); |
---|
211 | printf("nodeinfo: node %i\n", nodeInfo->nodeNumber()); |
---|
212 | { |
---|
213 | const CbcIntegerBranchingObject* brPrint = |
---|
214 | dynamic_cast<const CbcIntegerBranchingObject*>(nodeInfo->parentBranch()); |
---|
215 | if (!brPrint) { |
---|
216 | printf(" parentBranch: NULL\n"); |
---|
217 | } else { |
---|
218 | const double* downBounds = brPrint->downBounds(); |
---|
219 | const double* upBounds = brPrint->upBounds(); |
---|
220 | int variable = brPrint->variable(); |
---|
221 | int way = brPrint->way(); |
---|
222 | printf(" parentBranch: var %i downBd [%i,%i] upBd [%i,%i] way %i\n", |
---|
223 | variable, static_cast<int>(downBounds[0]), static_cast<int>(downBounds[1]), |
---|
224 | static_cast<int>(upBounds[0]), static_cast<int>(upBounds[1]), way); |
---|
225 | } |
---|
226 | } |
---|
227 | if (! node) { |
---|
228 | printf(" owner: NULL\n"); |
---|
229 | } else { |
---|
230 | printf(" owner: node %i depth %i onTree %i active %i", |
---|
231 | node->nodeNumber(), node->depth(), node->onTree(), node->active()); |
---|
232 | const OsiBranchingObject* osibr = |
---|
233 | nodeInfo->owner()->branchingObject(); |
---|
234 | const CbcBranchingObject* cbcbr = |
---|
235 | dynamic_cast<const CbcBranchingObject*>(osibr); |
---|
236 | const CbcIntegerBranchingObject* brPrint = |
---|
237 | dynamic_cast<const CbcIntegerBranchingObject*>(cbcbr); |
---|
238 | if (!brPrint) { |
---|
239 | printf(" ownerBranch: NULL\n"); |
---|
240 | } else { |
---|
241 | const double* downBounds = brPrint->downBounds(); |
---|
242 | const double* upBounds = brPrint->upBounds(); |
---|
243 | int variable = brPrint->variable(); |
---|
244 | int way = brPrint->way(); |
---|
245 | printf(" ownerbranch: var %i downBd [%i,%i] upBd [%i,%i] way %i\n", |
---|
246 | variable, static_cast<int>(downBounds[0]), static_cast<int>(downBounds[1]), |
---|
247 | static_cast<int>(upBounds[0]), static_cast<int>(upBounds[1]), way); |
---|
248 | } |
---|
249 | } |
---|
250 | nodeInfo = nodeInfo->parent(); |
---|
251 | } |
---|
252 | } |
---|
253 | |
---|
254 | void |
---|
255 | CbcHeuristic::debugNodes() |
---|
256 | { |
---|
257 | CbcHeurDebugNodes(model_); |
---|
258 | } |
---|
259 | |
---|
260 | void |
---|
261 | CbcHeuristic::printDistanceToNodes() |
---|
262 | { |
---|
263 | const CbcNode* currentNode = model_->currentNode(); |
---|
264 | if (currentNode != NULL) { |
---|
265 | CbcHeuristicNode* nodeDesc = new CbcHeuristicNode(*model_); |
---|
266 | for (int i = runNodes_.size() - 1; i >= 0; --i) { |
---|
267 | nodeDesc->distance(runNodes_.node(i)); |
---|
268 | } |
---|
269 | runNodes_.append(nodeDesc); |
---|
270 | } |
---|
271 | } |
---|
272 | |
---|
273 | bool |
---|
274 | CbcHeuristic::shouldHeurRun(int whereFrom) |
---|
275 | { |
---|
276 | assert (whereFrom >= 0 && whereFrom < 16); |
---|
277 | // take off 8 (code - likes new solution) |
---|
278 | whereFrom &= 7; |
---|
279 | if ((whereFrom_&(1 << whereFrom)) == 0) |
---|
280 | return false; |
---|
281 | // No longer used for original purpose - so use for ever run at all JJF |
---|
282 | #ifndef JJF_ONE |
---|
283 | // Don't run if hot start or no rows! |
---|
284 | if (model_ && (model_->hotstartSolution()||!model_->getNumRows())) |
---|
285 | return false; |
---|
286 | else |
---|
287 | return true; |
---|
288 | #else |
---|
289 | #ifdef JJF_ZERO |
---|
290 | const CbcNode* currentNode = model_->currentNode(); |
---|
291 | if (currentNode == NULL) { |
---|
292 | return false; |
---|
293 | } |
---|
294 | |
---|
295 | debugNodes(); |
---|
296 | // return false; |
---|
297 | |
---|
298 | const int depth = currentNode->depth(); |
---|
299 | #else |
---|
300 | int depth = model_->currentDepth(); |
---|
301 | #endif |
---|
302 | |
---|
303 | const int nodeCount = model_->getNodeCount(); // FIXME: check that this is |
---|
304 | // correct in parallel |
---|
305 | |
---|
306 | if (nodeCount == 0 || depth <= shallowDepth_) { |
---|
307 | // what to do when we are in the shallow part of the tree |
---|
308 | if (model_->getCurrentPassNumber() == 1) { |
---|
309 | // first time in the node... |
---|
310 | numInvocationsInShallow_ = 0; |
---|
311 | } |
---|
312 | ++numInvocationsInShallow_; |
---|
313 | // Very large howOftenShallow_ will give the original test: |
---|
314 | // (model_->getCurrentPassNumber() != 1) |
---|
315 | // if ((numInvocationsInShallow_ % howOftenShallow_) != 1) { |
---|
316 | if ((numInvocationsInShallow_ % howOftenShallow_) != 0) { |
---|
317 | return false; |
---|
318 | } |
---|
319 | // LL: should we save these nodes in the list of nodes where the heur was |
---|
320 | // LL: run? |
---|
321 | #ifndef JJF_ONE |
---|
322 | if (currentNode != NULL) { |
---|
323 | // Get where we are and create the appropriate CbcHeuristicNode object |
---|
324 | CbcHeuristicNode* nodeDesc = new CbcHeuristicNode(*model_); |
---|
325 | runNodes_.append(nodeDesc); |
---|
326 | } |
---|
327 | #endif |
---|
328 | } else { |
---|
329 | // deeper in the tree |
---|
330 | if (model_->getCurrentPassNumber() == 1) { |
---|
331 | // first time in the node... |
---|
332 | ++numInvocationsInDeep_; |
---|
333 | } |
---|
334 | if (numInvocationsInDeep_ - lastRunDeep_ < howOften_) { |
---|
335 | return false; |
---|
336 | } |
---|
337 | if (model_->getCurrentPassNumber() > 1) { |
---|
338 | // Run the heuristic only when first entering the node. |
---|
339 | // LL: I don't think this is right. It should run just before strong |
---|
340 | // LL: branching, I believe. |
---|
341 | return false; |
---|
342 | } |
---|
343 | // Get where we are and create the appropriate CbcHeuristicNode object |
---|
344 | CbcHeuristicNode* nodeDesc = new CbcHeuristicNode(*model_); |
---|
345 | //#ifdef PRINT_DEBUG |
---|
346 | #ifndef JJF_ONE |
---|
347 | const double minDistanceToRun = 1.5 * log((double)depth) / log((double)2); |
---|
348 | #else |
---|
349 | const double minDistanceToRun = minDistanceToRun_; |
---|
350 | #endif |
---|
351 | #ifdef PRINT_DEBUG |
---|
352 | double minDistance = nodeDesc->minDistance(runNodes_); |
---|
353 | std::cout << "minDistance = " << minDistance |
---|
354 | << ", minDistanceToRun = " << minDistanceToRun << std::endl; |
---|
355 | #endif |
---|
356 | if (nodeDesc->minDistanceIsSmall(runNodes_, minDistanceToRun)) { |
---|
357 | delete nodeDesc; |
---|
358 | return false; |
---|
359 | } |
---|
360 | runNodes_.append(nodeDesc); |
---|
361 | lastRunDeep_ = numInvocationsInDeep_; |
---|
362 | // ++lastRunDeep_; |
---|
363 | } |
---|
364 | ++numRuns_; |
---|
365 | return true; |
---|
366 | #endif |
---|
367 | } |
---|
368 | |
---|
369 | bool |
---|
370 | CbcHeuristic::shouldHeurRun_randomChoice() |
---|
371 | { |
---|
372 | if (!when_) |
---|
373 | return false; |
---|
374 | int depth = model_->currentDepth(); |
---|
375 | // when_ -999 is special marker to force to run |
---|
376 | if (depth != 0 && when_ != -999) { |
---|
377 | const double numerator = depth * depth; |
---|
378 | const double denominator = exp(depth * log(2.0)); |
---|
379 | double probability = numerator / denominator; |
---|
380 | double randomNumber = randomNumberGenerator_.randomDouble(); |
---|
381 | int when = when_ % 100; |
---|
382 | if (when > 2 && when < 8) { |
---|
383 | /* JJF adjustments |
---|
384 | 3 only at root and if no solution |
---|
385 | 4 only at root and if this heuristic has not got solution |
---|
386 | 5 decay (but only if no solution) |
---|
387 | 6 if depth <3 or decay |
---|
388 | 7 run up to 2 times if solution found 4 otherwise |
---|
389 | */ |
---|
390 | switch (when) { |
---|
391 | case 3: |
---|
392 | default: |
---|
393 | if (model_->bestSolution()) |
---|
394 | probability = -1.0; |
---|
395 | break; |
---|
396 | case 4: |
---|
397 | if (numberSolutionsFound_) |
---|
398 | probability = -1.0; |
---|
399 | break; |
---|
400 | case 5: |
---|
401 | assert (decayFactor_); |
---|
402 | if (model_->bestSolution()) { |
---|
403 | probability = -1.0; |
---|
404 | } else if (numCouldRun_ > 1000) { |
---|
405 | decayFactor_ *= 0.99; |
---|
406 | probability *= decayFactor_; |
---|
407 | } |
---|
408 | break; |
---|
409 | case 6: |
---|
410 | if (depth >= 3) { |
---|
411 | if ((numCouldRun_ % howOften_) == 0 && |
---|
412 | numberSolutionsFound_*howOften_ < numCouldRun_) { |
---|
413 | //#define COIN_DEVELOP |
---|
414 | #ifdef COIN_DEVELOP |
---|
415 | int old = howOften_; |
---|
416 | #endif |
---|
417 | howOften_ = CoinMin(CoinMax(static_cast<int> (howOften_ * 1.1), howOften_ + 1), 1000000); |
---|
418 | #ifdef COIN_DEVELOP |
---|
419 | printf("Howoften changed from %d to %d for %s\n", |
---|
420 | old, howOften_, heuristicName_.c_str()); |
---|
421 | #endif |
---|
422 | } |
---|
423 | probability = 1.0 / howOften_; |
---|
424 | if (model_->bestSolution()) |
---|
425 | probability *= 0.5; |
---|
426 | } else { |
---|
427 | probability = 1.1; |
---|
428 | } |
---|
429 | break; |
---|
430 | case 7: |
---|
431 | if ((model_->bestSolution() && numRuns_ >= 2) || numRuns_ >= 4) |
---|
432 | probability = -1.0; |
---|
433 | break; |
---|
434 | } |
---|
435 | } |
---|
436 | if (randomNumber > probability) |
---|
437 | return false; |
---|
438 | |
---|
439 | if (model_->getCurrentPassNumber() > 1) |
---|
440 | return false; |
---|
441 | #ifdef COIN_DEVELOP |
---|
442 | printf("Running %s, random %g probability %g\n", |
---|
443 | heuristicName_.c_str(), randomNumber, probability); |
---|
444 | #endif |
---|
445 | } else { |
---|
446 | #ifdef COIN_DEVELOP |
---|
447 | printf("Running %s, depth %d when %d\n", |
---|
448 | heuristicName_.c_str(), depth, when_); |
---|
449 | #endif |
---|
450 | } |
---|
451 | ++numRuns_; |
---|
452 | return true; |
---|
453 | } |
---|
454 | |
---|
455 | // Resets stuff if model changes |
---|
456 | void |
---|
457 | CbcHeuristic::resetModel(CbcModel * model) |
---|
458 | { |
---|
459 | model_ = model; |
---|
460 | } |
---|
461 | // Set seed |
---|
462 | void |
---|
463 | CbcHeuristic::setSeed(int value) |
---|
464 | { |
---|
465 | if (value==0) { |
---|
466 | double time = fabs(CoinGetTimeOfDay()); |
---|
467 | while (time>=COIN_INT_MAX) |
---|
468 | time *= 0.5; |
---|
469 | value = static_cast<int>(time); |
---|
470 | char printArray[100]; |
---|
471 | sprintf(printArray, "using time of day seed was changed from %d to %d", |
---|
472 | randomNumberGenerator_.getSeed(), value); |
---|
473 | if (model_) |
---|
474 | model_->messageHandler()->message(CBC_FPUMP1, model_->messages()) |
---|
475 | << printArray |
---|
476 | << CoinMessageEol; |
---|
477 | } |
---|
478 | randomNumberGenerator_.setSeed(value); |
---|
479 | } |
---|
480 | // Get seed |
---|
481 | int |
---|
482 | CbcHeuristic::getSeed() const |
---|
483 | { |
---|
484 | return randomNumberGenerator_.getSeed(); |
---|
485 | } |
---|
486 | |
---|
487 | // Create C++ lines to get to current state |
---|
488 | void |
---|
489 | CbcHeuristic::generateCpp( FILE * fp, const char * heuristic) |
---|
490 | { |
---|
491 | // hard coded as CbcHeuristic virtual |
---|
492 | if (when_ != 2) |
---|
493 | fprintf(fp, "3 %s.setWhen(%d);\n", heuristic, when_); |
---|
494 | else |
---|
495 | fprintf(fp, "4 %s.setWhen(%d);\n", heuristic, when_); |
---|
496 | if (numberNodes_ != 200) |
---|
497 | fprintf(fp, "3 %s.setNumberNodes(%d);\n", heuristic, numberNodes_); |
---|
498 | else |
---|
499 | fprintf(fp, "4 %s.setNumberNodes(%d);\n", heuristic, numberNodes_); |
---|
500 | if (feasibilityPumpOptions_ != -1) |
---|
501 | fprintf(fp, "3 %s.setFeasibilityPumpOptions(%d);\n", heuristic, feasibilityPumpOptions_); |
---|
502 | else |
---|
503 | fprintf(fp, "4 %s.setFeasibilityPumpOptions(%d);\n", heuristic, feasibilityPumpOptions_); |
---|
504 | if (fractionSmall_ != 1.0) |
---|
505 | fprintf(fp, "3 %s.setFractionSmall(%g);\n", heuristic, fractionSmall_); |
---|
506 | else |
---|
507 | fprintf(fp, "4 %s.setFractionSmall(%g);\n", heuristic, fractionSmall_); |
---|
508 | if (heuristicName_ != "Unknown") |
---|
509 | fprintf(fp, "3 %s.setHeuristicName(\"%s\");\n", |
---|
510 | heuristic, heuristicName_.c_str()) ; |
---|
511 | else |
---|
512 | fprintf(fp, "4 %s.setHeuristicName(\"%s\");\n", |
---|
513 | heuristic, heuristicName_.c_str()) ; |
---|
514 | if (decayFactor_ != 0.0) |
---|
515 | fprintf(fp, "3 %s.setDecayFactor(%g);\n", heuristic, decayFactor_); |
---|
516 | else |
---|
517 | fprintf(fp, "4 %s.setDecayFactor(%g);\n", heuristic, decayFactor_); |
---|
518 | if (switches_ != 0) |
---|
519 | fprintf(fp, "3 %s.setSwitches(%d);\n", heuristic, switches_); |
---|
520 | else |
---|
521 | fprintf(fp, "4 %s.setSwitches(%d);\n", heuristic, switches_); |
---|
522 | if (whereFrom_ != DEFAULT_WHERE) |
---|
523 | fprintf(fp, "3 %s.setWhereFrom(%d);\n", heuristic, whereFrom_); |
---|
524 | else |
---|
525 | fprintf(fp, "4 %s.setWhereFrom(%d);\n", heuristic, whereFrom_); |
---|
526 | if (shallowDepth_ != 1) |
---|
527 | fprintf(fp, "3 %s.setShallowDepth(%d);\n", heuristic, shallowDepth_); |
---|
528 | else |
---|
529 | fprintf(fp, "4 %s.setShallowDepth(%d);\n", heuristic, shallowDepth_); |
---|
530 | if (howOftenShallow_ != 1) |
---|
531 | fprintf(fp, "3 %s.setHowOftenShallow(%d);\n", heuristic, howOftenShallow_); |
---|
532 | else |
---|
533 | fprintf(fp, "4 %s.setHowOftenShallow(%d);\n", heuristic, howOftenShallow_); |
---|
534 | if (minDistanceToRun_ != 1) |
---|
535 | fprintf(fp, "3 %s.setMinDistanceToRun(%d);\n", heuristic, minDistanceToRun_); |
---|
536 | else |
---|
537 | fprintf(fp, "4 %s.setMinDistanceToRun(%d);\n", heuristic, minDistanceToRun_); |
---|
538 | } |
---|
539 | // Destructor |
---|
540 | CbcHeuristic::~CbcHeuristic () |
---|
541 | { |
---|
542 | delete [] inputSolution_; |
---|
543 | } |
---|
544 | |
---|
545 | // update model |
---|
546 | void CbcHeuristic::setModel(CbcModel * model) |
---|
547 | { |
---|
548 | model_ = model; |
---|
549 | } |
---|
550 | /* Clone but .. |
---|
551 | type 0 clone solver, 1 clone continuous solver |
---|
552 | Add 2 to say without integer variables which are at low priority |
---|
553 | Add 4 to say quite likely infeasible so give up easily.*/ |
---|
554 | OsiSolverInterface * |
---|
555 | CbcHeuristic::cloneBut(int type) |
---|
556 | { |
---|
557 | OsiSolverInterface * solver; |
---|
558 | if ((type&1) == 0 || !model_->continuousSolver()) |
---|
559 | solver = model_->solver()->clone(); |
---|
560 | else |
---|
561 | solver = model_->continuousSolver()->clone(); |
---|
562 | #ifdef COIN_HAS_CLP |
---|
563 | OsiClpSolverInterface * clpSolver |
---|
564 | = dynamic_cast<OsiClpSolverInterface *> (solver); |
---|
565 | #endif |
---|
566 | if ((type&2) != 0) { |
---|
567 | int n = model_->numberObjects(); |
---|
568 | int priority = model_->continuousPriority(); |
---|
569 | if (priority < COIN_INT_MAX) { |
---|
570 | for (int i = 0; i < n; i++) { |
---|
571 | const OsiObject * obj = model_->object(i); |
---|
572 | const CbcSimpleInteger * thisOne = |
---|
573 | dynamic_cast <const CbcSimpleInteger *> (obj); |
---|
574 | if (thisOne) { |
---|
575 | int iColumn = thisOne->columnNumber(); |
---|
576 | if (thisOne->priority() >= priority) |
---|
577 | solver->setContinuous(iColumn); |
---|
578 | } |
---|
579 | } |
---|
580 | } |
---|
581 | #ifdef COIN_HAS_CLP |
---|
582 | if (clpSolver) { |
---|
583 | for (int i = 0; i < n; i++) { |
---|
584 | const OsiObject * obj = model_->object(i); |
---|
585 | const CbcSimpleInteger * thisOne = |
---|
586 | dynamic_cast <const CbcSimpleInteger *> (obj); |
---|
587 | if (thisOne) { |
---|
588 | int iColumn = thisOne->columnNumber(); |
---|
589 | if (clpSolver->isOptionalInteger(iColumn)) |
---|
590 | clpSolver->setContinuous(iColumn); |
---|
591 | } |
---|
592 | } |
---|
593 | } |
---|
594 | #endif |
---|
595 | } |
---|
596 | #ifdef COIN_HAS_CLP |
---|
597 | if ((type&4) != 0 && clpSolver) { |
---|
598 | int options = clpSolver->getModelPtr()->moreSpecialOptions(); |
---|
599 | clpSolver->getModelPtr()->setMoreSpecialOptions(options | 64); |
---|
600 | } |
---|
601 | #endif |
---|
602 | return solver; |
---|
603 | } |
---|
604 | // Whether to exit at once on gap |
---|
605 | bool |
---|
606 | CbcHeuristic::exitNow(double bestObjective) const |
---|
607 | { |
---|
608 | if ((switches_&2048) != 0) { |
---|
609 | // exit may be forced - but unset for next time |
---|
610 | switches_ &= ~2048; |
---|
611 | if ((switches_&1024) != 0) |
---|
612 | return true; |
---|
613 | } else if ((switches_&1) == 0) { |
---|
614 | return false; |
---|
615 | } |
---|
616 | // See if can stop on gap |
---|
617 | OsiSolverInterface * solver = model_->solver(); |
---|
618 | double bestPossibleObjective = solver->getObjValue() * solver->getObjSense(); |
---|
619 | double absGap = CoinMax(model_->getAllowableGap(), |
---|
620 | model_->getHeuristicGap()); |
---|
621 | double fracGap = CoinMax(model_->getAllowableFractionGap(), |
---|
622 | model_->getHeuristicFractionGap()); |
---|
623 | double testGap = CoinMax(absGap, fracGap * |
---|
624 | CoinMax(fabs(bestObjective), |
---|
625 | fabs(bestPossibleObjective))); |
---|
626 | |
---|
627 | if (bestObjective - bestPossibleObjective < testGap |
---|
628 | && model_->getCutoffIncrement() >= 0.0) { |
---|
629 | return true; |
---|
630 | } else { |
---|
631 | return false; |
---|
632 | } |
---|
633 | } |
---|
634 | #ifdef HISTORY_STATISTICS |
---|
635 | extern bool getHistoryStatistics_; |
---|
636 | #endif |
---|
637 | static double sizeRatio(int numberRowsNow, int numberColumnsNow, |
---|
638 | int numberRowsStart, int numberColumnsStart) |
---|
639 | { |
---|
640 | double valueNow; |
---|
641 | if (numberRowsNow*10 > numberColumnsNow || numberColumnsNow < 200) { |
---|
642 | valueNow = 2 * numberRowsNow + numberColumnsNow; |
---|
643 | } else { |
---|
644 | // long and thin - rows are more important |
---|
645 | if (numberRowsNow*40 > numberColumnsNow) |
---|
646 | valueNow = 10 * numberRowsNow + numberColumnsNow; |
---|
647 | else |
---|
648 | valueNow = 200 * numberRowsNow + numberColumnsNow; |
---|
649 | } |
---|
650 | double valueStart; |
---|
651 | if (numberRowsStart*10 > numberColumnsStart || numberColumnsStart < 200) { |
---|
652 | valueStart = 2 * numberRowsStart + numberColumnsStart; |
---|
653 | } else { |
---|
654 | // long and thin - rows are more important |
---|
655 | if (numberRowsStart*40 > numberColumnsStart) |
---|
656 | valueStart = 10 * numberRowsStart + numberColumnsStart; |
---|
657 | else |
---|
658 | valueStart = 200 * numberRowsStart + numberColumnsStart; |
---|
659 | } |
---|
660 | //printf("sizeProblem Now %g, %d rows, %d columns\nsizeProblem Start %g, %d rows, %d columns\n", |
---|
661 | // valueNow,numberRowsNow,numberColumnsNow, |
---|
662 | // valueStart,numberRowsStart,numberColumnsStart); |
---|
663 | if (10*numberRowsNow < 8*numberRowsStart || 10*numberColumnsNow < 7*numberColumnsStart) |
---|
664 | return valueNow / valueStart; |
---|
665 | else if (10*numberRowsNow < 9*numberRowsStart) |
---|
666 | return 1.1*(valueNow / valueStart); |
---|
667 | else if (numberRowsNow < numberRowsStart) |
---|
668 | return 1.5*(valueNow / valueStart); |
---|
669 | else |
---|
670 | return 2.0*(valueNow / valueStart); |
---|
671 | } |
---|
672 | |
---|
673 | //static int saveModel=0; |
---|
674 | // Do mini branch and bound (return 1 if solution) |
---|
675 | int |
---|
676 | CbcHeuristic::smallBranchAndBound(OsiSolverInterface * solver, int numberNodes, |
---|
677 | double * newSolution, double & newSolutionValue, |
---|
678 | double cutoff, std::string name) const |
---|
679 | { |
---|
680 | CbcEventHandler *eventHandler = model_->getEventHandler() ; |
---|
681 | // Use this fraction |
---|
682 | double fractionSmall = fractionSmall_; |
---|
683 | int maximumSolutions = model_->getMaximumSolutions(); |
---|
684 | int iterationMultiplier = 100; |
---|
685 | if (eventHandler) { |
---|
686 | typedef struct { |
---|
687 | double fractionSmall; |
---|
688 | double spareDouble[3]; |
---|
689 | OsiSolverInterface * solver; |
---|
690 | void * sparePointer[2]; |
---|
691 | int numberNodes; |
---|
692 | int maximumSolutions; |
---|
693 | int iterationMultiplier; |
---|
694 | int howOften; |
---|
695 | int spareInt[3]; |
---|
696 | } SmallMod; |
---|
697 | SmallMod temp; |
---|
698 | temp.solver=solver; |
---|
699 | temp.fractionSmall=fractionSmall; |
---|
700 | temp.numberNodes=numberNodes; |
---|
701 | temp.iterationMultiplier=iterationMultiplier; |
---|
702 | temp.howOften=howOften_; |
---|
703 | temp.maximumSolutions=maximumSolutions; |
---|
704 | CbcEventHandler::CbcAction status = |
---|
705 | eventHandler->event(CbcEventHandler::smallBranchAndBound, |
---|
706 | &temp); |
---|
707 | if (status==CbcEventHandler::killSolution) |
---|
708 | return -1; |
---|
709 | if (status==CbcEventHandler::takeAction) { |
---|
710 | fractionSmall=temp.fractionSmall; |
---|
711 | numberNodes=temp.numberNodes; |
---|
712 | iterationMultiplier=temp.iterationMultiplier; |
---|
713 | howOften_=temp.howOften; |
---|
714 | maximumSolutions=temp.maximumSolutions; |
---|
715 | } |
---|
716 | } |
---|
717 | #if 0 |
---|
718 | if (saveModel || model_->getMaximumSolutions()==100) { |
---|
719 | printf("writing model\n"); |
---|
720 | solver->writeMpsNative("before.mps", NULL, NULL, 2, 1); |
---|
721 | } |
---|
722 | #endif |
---|
723 | // size before |
---|
724 | int shiftRows = 0; |
---|
725 | if (numberNodes < 0) |
---|
726 | shiftRows = solver->getNumRows() - numberNodes_; |
---|
727 | int numberRowsStart = solver->getNumRows() - shiftRows; |
---|
728 | int numberColumnsStart = solver->getNumCols(); |
---|
729 | #ifdef CLP_INVESTIGATE |
---|
730 | printf("%s has %d rows, %d columns\n", |
---|
731 | name.c_str(), solver->getNumRows(), solver->getNumCols()); |
---|
732 | #endif |
---|
733 | double before = 2 * numberRowsStart + numberColumnsStart; |
---|
734 | if (before > 40000.0) { |
---|
735 | // fairly large - be more conservative |
---|
736 | double multiplier = 1.0 - 0.3 * CoinMin(100000.0, before - 40000.0) / 100000.0; |
---|
737 | if (multiplier < 1.0) { |
---|
738 | fractionSmall *= multiplier; |
---|
739 | #ifdef CLP_INVESTIGATE |
---|
740 | printf("changing fractionSmall from %g to %g for %s\n", |
---|
741 | fractionSmall_, fractionSmall, name.c_str()); |
---|
742 | #endif |
---|
743 | } |
---|
744 | } |
---|
745 | #ifdef COIN_HAS_CLP |
---|
746 | OsiClpSolverInterface * osiclp = dynamic_cast< OsiClpSolverInterface*> (solver); |
---|
747 | if (osiclp && (osiclp->specialOptions()&65536) == 0) { |
---|
748 | // go faster stripes |
---|
749 | if (osiclp->getNumRows() < 300 && osiclp->getNumCols() < 500) { |
---|
750 | osiclp->setupForRepeatedUse(2, 0); |
---|
751 | } else { |
---|
752 | osiclp->setupForRepeatedUse(0, 0); |
---|
753 | } |
---|
754 | // Turn this off if you get problems |
---|
755 | // Used to be automatically set |
---|
756 | osiclp->setSpecialOptions(osiclp->specialOptions() | (128 + 64 - 128)); |
---|
757 | ClpSimplex * lpSolver = osiclp->getModelPtr(); |
---|
758 | lpSolver->setSpecialOptions(lpSolver->specialOptions() | 0x01000000); // say is Cbc (and in branch and bound) |
---|
759 | lpSolver->setSpecialOptions(lpSolver->specialOptions() | |
---|
760 | (/*16384+*/4096 + 512 + 128)); |
---|
761 | } |
---|
762 | #endif |
---|
763 | #ifdef HISTORY_STATISTICS |
---|
764 | getHistoryStatistics_ = false; |
---|
765 | #endif |
---|
766 | #ifdef COIN_DEVELOP |
---|
767 | int status = 0; |
---|
768 | #endif |
---|
769 | int logLevel = model_->logLevel(); |
---|
770 | #define LEN_PRINT 250 |
---|
771 | char generalPrint[LEN_PRINT]; |
---|
772 | // Do presolve to see if possible |
---|
773 | int numberColumns = solver->getNumCols(); |
---|
774 | char * reset = NULL; |
---|
775 | int returnCode = 1; |
---|
776 | int saveModelOptions = model_->specialOptions(); |
---|
777 | //assert ((saveModelOptions&2048) == 0); |
---|
778 | model_->setSpecialOptions(saveModelOptions | 2048); |
---|
779 | if (fractionSmall<1.0) { |
---|
780 | int saveLogLevel = solver->messageHandler()->logLevel(); |
---|
781 | if (saveLogLevel == 1) |
---|
782 | solver->messageHandler()->setLogLevel(0); |
---|
783 | OsiPresolve * pinfo = new OsiPresolve(); |
---|
784 | int presolveActions = 0; |
---|
785 | // Allow dual stuff on integers |
---|
786 | presolveActions = 1; |
---|
787 | // Do not allow all +1 to be tampered with |
---|
788 | //if (allPlusOnes) |
---|
789 | //presolveActions |= 2; |
---|
790 | // allow transfer of costs |
---|
791 | // presolveActions |= 4; |
---|
792 | pinfo->setPresolveActions(presolveActions); |
---|
793 | OsiSolverInterface * presolvedModel = pinfo->presolvedModel(*solver, 1.0e-8, true, 2); |
---|
794 | delete pinfo; |
---|
795 | // see if too big |
---|
796 | |
---|
797 | if (presolvedModel) { |
---|
798 | int afterRows = presolvedModel->getNumRows(); |
---|
799 | int afterCols = presolvedModel->getNumCols(); |
---|
800 | //#define COIN_DEVELOP |
---|
801 | #ifdef COIN_DEVELOP_z |
---|
802 | if (numberNodes < 0) { |
---|
803 | solver->writeMpsNative("before.mps", NULL, NULL, 2, 1); |
---|
804 | presolvedModel->writeMpsNative("after1.mps", NULL, NULL, 2, 1); |
---|
805 | } |
---|
806 | #endif |
---|
807 | delete presolvedModel; |
---|
808 | double ratio = sizeRatio(afterRows - shiftRows, afterCols, |
---|
809 | numberRowsStart, numberColumnsStart); |
---|
810 | double after = 2 * afterRows + afterCols; |
---|
811 | if (ratio > fractionSmall && after > 300 && numberNodes >= 0) { |
---|
812 | // Need code to try again to compress further using used |
---|
813 | const int * used = model_->usedInSolution(); |
---|
814 | int maxUsed = 0; |
---|
815 | int iColumn; |
---|
816 | const double * lower = solver->getColLower(); |
---|
817 | const double * upper = solver->getColUpper(); |
---|
818 | for (iColumn = 0; iColumn < numberColumns; iColumn++) { |
---|
819 | if (upper[iColumn] > lower[iColumn]) { |
---|
820 | if (solver->isBinary(iColumn)) |
---|
821 | maxUsed = CoinMax(maxUsed, used[iColumn]); |
---|
822 | } |
---|
823 | } |
---|
824 | if (maxUsed) { |
---|
825 | reset = new char [numberColumns]; |
---|
826 | int nFix = 0; |
---|
827 | for (iColumn = 0; iColumn < numberColumns; iColumn++) { |
---|
828 | reset[iColumn] = 0; |
---|
829 | if (upper[iColumn] > lower[iColumn]) { |
---|
830 | if (solver->isBinary(iColumn) && used[iColumn] == maxUsed) { |
---|
831 | bool setValue = true; |
---|
832 | if (maxUsed == 1) { |
---|
833 | double randomNumber = randomNumberGenerator_.randomDouble(); |
---|
834 | if (randomNumber > 0.3) |
---|
835 | setValue = false; |
---|
836 | } |
---|
837 | if (setValue) { |
---|
838 | reset[iColumn] = 1; |
---|
839 | solver->setColLower(iColumn, 1.0); |
---|
840 | nFix++; |
---|
841 | } |
---|
842 | } |
---|
843 | } |
---|
844 | } |
---|
845 | pinfo = new OsiPresolve(); |
---|
846 | presolveActions = 0; |
---|
847 | // Allow dual stuff on integers |
---|
848 | presolveActions = 1; |
---|
849 | // Do not allow all +1 to be tampered with |
---|
850 | //if (allPlusOnes) |
---|
851 | //presolveActions |= 2; |
---|
852 | // allow transfer of costs |
---|
853 | // presolveActions |= 4; |
---|
854 | pinfo->setPresolveActions(presolveActions); |
---|
855 | presolvedModel = pinfo->presolvedModel(*solver, 1.0e-8, true, 2); |
---|
856 | delete pinfo; |
---|
857 | if (presolvedModel) { |
---|
858 | // see if too big |
---|
859 | int afterRows2 = presolvedModel->getNumRows(); |
---|
860 | int afterCols2 = presolvedModel->getNumCols(); |
---|
861 | delete presolvedModel; |
---|
862 | double ratio = sizeRatio(afterRows2 - shiftRows, afterCols2, |
---|
863 | numberRowsStart, numberColumnsStart); |
---|
864 | double after = 2 * afterRows2 + afterCols2; |
---|
865 | if (ratio > fractionSmall && (after > 300 || numberNodes < 0)) { |
---|
866 | sprintf(generalPrint, "Full problem %d rows %d columns, reduced to %d rows %d columns - %d fixed gives %d, %d - still too large", |
---|
867 | solver->getNumRows(), solver->getNumCols(), |
---|
868 | afterRows, afterCols, nFix, afterRows2, afterCols2); |
---|
869 | // If much too big - give up |
---|
870 | if (ratio > 0.75) |
---|
871 | returnCode = -1; |
---|
872 | } else { |
---|
873 | sprintf(generalPrint, "Full problem %d rows %d columns, reduced to %d rows %d columns - %d fixed gives %d, %d - ok now", |
---|
874 | solver->getNumRows(), solver->getNumCols(), |
---|
875 | afterRows, afterCols, nFix, afterRows2, afterCols2); |
---|
876 | } |
---|
877 | model_->messageHandler()->message(CBC_FPUMP1, model_->messages()) |
---|
878 | << generalPrint |
---|
879 | << CoinMessageEol; |
---|
880 | } else { |
---|
881 | returnCode = 2; // infeasible |
---|
882 | } |
---|
883 | } |
---|
884 | } else if (ratio > fractionSmall && after > 300 && numberNodes >=0) { |
---|
885 | returnCode = -1; |
---|
886 | } |
---|
887 | } else { |
---|
888 | returnCode = 2; // infeasible |
---|
889 | } |
---|
890 | solver->messageHandler()->setLogLevel(saveLogLevel); |
---|
891 | } |
---|
892 | if (returnCode == 2 || returnCode == -1) { |
---|
893 | model_->setSpecialOptions(saveModelOptions); |
---|
894 | delete [] reset; |
---|
895 | #ifdef HISTORY_STATISTICS |
---|
896 | getHistoryStatistics_ = true; |
---|
897 | #endif |
---|
898 | //printf("small no good\n"); |
---|
899 | return returnCode; |
---|
900 | } |
---|
901 | // Reduce printout |
---|
902 | bool takeHint; |
---|
903 | OsiHintStrength strength; |
---|
904 | solver->getHintParam(OsiDoReducePrint, takeHint, strength); |
---|
905 | solver->setHintParam(OsiDoReducePrint, true, OsiHintTry); |
---|
906 | solver->setHintParam(OsiDoPresolveInInitial, false, OsiHintTry); |
---|
907 | double signedCutoff = cutoff*solver->getObjSense(); |
---|
908 | solver->setDblParam(OsiDualObjectiveLimit, signedCutoff); |
---|
909 | solver->initialSolve(); |
---|
910 | if (solver->isProvenOptimal()) { |
---|
911 | CglPreProcess process; |
---|
912 | OsiSolverInterface * solver2 = NULL; |
---|
913 | if ((model_->moreSpecialOptions()&65536)!=0) |
---|
914 | process.setOptions(2+4+8+16); // no cuts |
---|
915 | else |
---|
916 | process.setOptions(16); // no complicated dupcol stuff |
---|
917 | /* Do not try and produce equality cliques and |
---|
918 | do up to 2 passes (normally) 5 if restart */ |
---|
919 | int numberPasses = 2; |
---|
920 | if ((model_->moreSpecialOptions2()&16)!=0) { |
---|
921 | // quick |
---|
922 | process.setOptions(2+4+8+16); // no cuts |
---|
923 | numberPasses = 1; |
---|
924 | } |
---|
925 | if (numberNodes < 0) { |
---|
926 | numberPasses = 5; |
---|
927 | // Say some rows cuts |
---|
928 | int numberRows = solver->getNumRows(); |
---|
929 | if (numberNodes_ < numberRows && true /* think */) { |
---|
930 | char * type = new char[numberRows]; |
---|
931 | memset(type, 0, numberNodes_); |
---|
932 | memset(type + numberNodes_, 1, numberRows - numberNodes_); |
---|
933 | process.passInRowTypes(type, numberRows); |
---|
934 | delete [] type; |
---|
935 | } |
---|
936 | } |
---|
937 | if (logLevel <= 1) |
---|
938 | process.messageHandler()->setLogLevel(0); |
---|
939 | if (!solver->defaultHandler()&& |
---|
940 | solver->messageHandler()->logLevel(0)!=-1000) |
---|
941 | process.passInMessageHandler(solver->messageHandler()); |
---|
942 | #ifdef CGL_DEBUG |
---|
943 | /* |
---|
944 | We're debugging. (specialOptions 1) |
---|
945 | */ |
---|
946 | if ((model_->specialOptions()&1) != 0) { |
---|
947 | const OsiRowCutDebugger *debugger = solver->getRowCutDebugger() ; |
---|
948 | if (debugger) { |
---|
949 | process.setApplicationData(const_cast<double *>(debugger->optimalSolution())); |
---|
950 | } |
---|
951 | } |
---|
952 | #endif |
---|
953 | solver2 = process.preProcessNonDefault(*solver, false, |
---|
954 | numberPasses); |
---|
955 | if (!solver2) { |
---|
956 | if (logLevel > 1) |
---|
957 | printf("Pre-processing says infeasible\n"); |
---|
958 | returnCode = 2; // so will be infeasible |
---|
959 | } else { |
---|
960 | #ifdef COIN_DEVELOP_z |
---|
961 | if (numberNodes < 0) { |
---|
962 | solver2->writeMpsNative("after2.mps", NULL, NULL, 2, 1); |
---|
963 | } |
---|
964 | #endif |
---|
965 | // see if too big |
---|
966 | double ratio = sizeRatio(solver2->getNumRows() - shiftRows, solver2->getNumCols(), |
---|
967 | numberRowsStart, numberColumnsStart); |
---|
968 | double after = 2 * solver2->getNumRows() + solver2->getNumCols(); |
---|
969 | if (ratio > fractionSmall && (after > 300 || numberNodes < 0)) { |
---|
970 | sprintf(generalPrint, "Full problem %d rows %d columns, reduced to %d rows %d columns - too large", |
---|
971 | solver->getNumRows(), solver->getNumCols(), |
---|
972 | solver2->getNumRows(), solver2->getNumCols()); |
---|
973 | model_->messageHandler()->message(CBC_FPUMP1, model_->messages()) |
---|
974 | << generalPrint |
---|
975 | << CoinMessageEol; |
---|
976 | returnCode = -1; |
---|
977 | //printf("small no good2\n"); |
---|
978 | } else { |
---|
979 | sprintf(generalPrint, "Full problem %d rows %d columns, reduced to %d rows %d columns", |
---|
980 | solver->getNumRows(), solver->getNumCols(), |
---|
981 | solver2->getNumRows(), solver2->getNumCols()); |
---|
982 | model_->messageHandler()->message(CBC_FPUMP1, model_->messages()) |
---|
983 | << generalPrint |
---|
984 | << CoinMessageEol; |
---|
985 | } |
---|
986 | #ifdef CGL_DEBUG |
---|
987 | if ((model_->specialOptions()&1) != 0) { |
---|
988 | const OsiRowCutDebugger *debugger = solver2->getRowCutDebugger() ; |
---|
989 | if (debugger) { |
---|
990 | printf("On optimal path after preprocessing\n"); |
---|
991 | } |
---|
992 | } |
---|
993 | #endif |
---|
994 | if (returnCode == 1) { |
---|
995 | solver2->resolve(); |
---|
996 | CbcModel model(*solver2); |
---|
997 | double startTime=model_->getDblParam(CbcModel::CbcStartSeconds); |
---|
998 | model.setDblParam(CbcModel::CbcStartSeconds,startTime); |
---|
999 | // move seed across |
---|
1000 | model.randomNumberGenerator()->setSeed(model_->randomNumberGenerator()->getSeed()); |
---|
1001 | if (numberNodes >= 0) { |
---|
1002 | // normal |
---|
1003 | model.setSpecialOptions(saveModelOptions | 2048); |
---|
1004 | if (logLevel <= 1 && feasibilityPumpOptions_ != -3) |
---|
1005 | model.setLogLevel(0); |
---|
1006 | else |
---|
1007 | model.setLogLevel(logLevel); |
---|
1008 | // No small fathoming |
---|
1009 | model.setFastNodeDepth(-1); |
---|
1010 | model.setCutoff(signedCutoff); |
---|
1011 | model.setStrongStrategy(0); |
---|
1012 | // Don't do if original fraction > 1.0 and too large |
---|
1013 | if (fractionSmall_>1.0 && fractionSmall_ < 1000000.0) { |
---|
1014 | /* 1.4 means -1 nodes if >.4 |
---|
1015 | 2.4 means -1 nodes if >.5 and 0 otherwise |
---|
1016 | 3.4 means -1 nodes if >.6 and 0 or 5 |
---|
1017 | 4.4 means -1 nodes if >.7 and 0, 5 or 10 |
---|
1018 | */ |
---|
1019 | double fraction = fractionSmall_-floor(fractionSmall_); |
---|
1020 | if (ratio>fraction) { |
---|
1021 | int type = static_cast<int>(floor(fractionSmall_*0.1)); |
---|
1022 | int over = static_cast<int>(ceil(ratio-fraction)); |
---|
1023 | int maxNodes[]={-1,0,5,10}; |
---|
1024 | if (type>over) |
---|
1025 | numberNodes=maxNodes[type-over]; |
---|
1026 | else |
---|
1027 | numberNodes=-1; |
---|
1028 | } |
---|
1029 | } |
---|
1030 | model.setMaximumNodes(numberNodes); |
---|
1031 | model.solver()->setHintParam(OsiDoReducePrint, true, OsiHintTry); |
---|
1032 | if ((saveModelOptions&2048) == 0) |
---|
1033 | model.setMoreSpecialOptions(model_->moreSpecialOptions()); |
---|
1034 | model.setMoreSpecialOptions2(model_->moreSpecialOptions2()); |
---|
1035 | // off conflict analysis |
---|
1036 | model.setMoreSpecialOptions(model.moreSpecialOptions()&~4194304); |
---|
1037 | |
---|
1038 | // Lightweight |
---|
1039 | CbcStrategyDefaultSubTree strategy(model_, 1, 5, 1, 0); |
---|
1040 | model.setStrategy(strategy); |
---|
1041 | model.solver()->setIntParam(OsiMaxNumIterationHotStart, 10); |
---|
1042 | model.setMaximumCutPassesAtRoot(CoinMin(20, CoinAbs(model_->getMaximumCutPassesAtRoot()))); |
---|
1043 | model.setMaximumCutPasses(CoinMin(10, model_->getMaximumCutPasses())); |
---|
1044 | // Set best solution (even if bad for this submodel) |
---|
1045 | if (model_->bestSolution()) { |
---|
1046 | const double * bestSolution = model_->bestSolution(); |
---|
1047 | int numberColumns2 = model.solver()->getNumCols(); |
---|
1048 | double * bestSolution2 = new double [numberColumns2]; |
---|
1049 | const int * originalColumns = process.originalColumns(); |
---|
1050 | for (int iColumn=0;iColumn<numberColumns2;iColumn++) { |
---|
1051 | int jColumn = originalColumns[iColumn]; |
---|
1052 | bestSolution2[iColumn] = bestSolution[jColumn]; |
---|
1053 | } |
---|
1054 | model.setBestSolution(bestSolution2,numberColumns2, |
---|
1055 | 1.0e50, |
---|
1056 | false); |
---|
1057 | model.setSolutionCount(1); |
---|
1058 | maximumSolutions++; |
---|
1059 | delete [] bestSolution2; |
---|
1060 | } |
---|
1061 | } else { |
---|
1062 | model.setSpecialOptions(saveModelOptions); |
---|
1063 | model_->messageHandler()->message(CBC_RESTART, model_->messages()) |
---|
1064 | << solver2->getNumRows() << solver2->getNumCols() |
---|
1065 | << CoinMessageEol; |
---|
1066 | // going for full search and copy across more stuff |
---|
1067 | model.gutsOfCopy(*model_, 2); |
---|
1068 | #ifdef CGL_DEBUG |
---|
1069 | if ((model_->specialOptions()&1) != 0) { |
---|
1070 | const OsiRowCutDebugger *debugger = model.solver()->getRowCutDebugger() ; |
---|
1071 | if (debugger) { |
---|
1072 | printf("On optimal path BB\n"); |
---|
1073 | } |
---|
1074 | } |
---|
1075 | #endif |
---|
1076 | assert (!model_->heuristicModel()); |
---|
1077 | model_->setHeuristicModel(&model); |
---|
1078 | for (int i = 0; i < model.numberCutGenerators(); i++) { |
---|
1079 | CbcCutGenerator * generator = model.cutGenerator(i); |
---|
1080 | CglGomory * gomory = dynamic_cast<CglGomory *> |
---|
1081 | (generator->generator()); |
---|
1082 | if (gomory&&gomory->originalSolver()) |
---|
1083 | gomory->passInOriginalSolver(model.solver()); |
---|
1084 | generator->setTiming(true); |
---|
1085 | // Turn on if was turned on |
---|
1086 | int iOften = model_->cutGenerator(i)->howOften(); |
---|
1087 | #ifdef CLP_INVESTIGATE |
---|
1088 | printf("Gen %d often %d %d\n", |
---|
1089 | i, generator->howOften(), |
---|
1090 | iOften); |
---|
1091 | #endif |
---|
1092 | if (iOften > 0) |
---|
1093 | generator->setHowOften(iOften % 1000000); |
---|
1094 | if (model_->cutGenerator(i)->howOftenInSub() == -200) |
---|
1095 | generator->setHowOften(-100); |
---|
1096 | } |
---|
1097 | model.setCutoff(signedCutoff); |
---|
1098 | // make sure can't do nested search! but allow heuristics |
---|
1099 | model.setSpecialOptions((model.specialOptions()&(~(512 + 2048))) | 1024); |
---|
1100 | bool takeHint; |
---|
1101 | OsiHintStrength strength; |
---|
1102 | // Switch off printing if asked to |
---|
1103 | model_->solver()->getHintParam(OsiDoReducePrint, takeHint, strength); |
---|
1104 | model.solver()->setHintParam(OsiDoReducePrint, takeHint, strength); |
---|
1105 | // no cut generators if none in parent |
---|
1106 | CbcStrategyDefault |
---|
1107 | strategy(model_->numberCutGenerators() ? 1 : -1, |
---|
1108 | model_->numberStrong(), |
---|
1109 | model_->numberBeforeTrust()); |
---|
1110 | // Set up pre-processing - no |
---|
1111 | strategy.setupPreProcessing(0); // was (4); |
---|
1112 | model.setStrategy(strategy); |
---|
1113 | //model.solver()->writeMps("crunched"); |
---|
1114 | int numberCuts = process.cuts().sizeRowCuts(); |
---|
1115 | if (numberCuts) { |
---|
1116 | // add in cuts |
---|
1117 | CglStored cuts = process.cuts(); |
---|
1118 | model.addCutGenerator(&cuts, 1, "Stored from first"); |
---|
1119 | model.cutGenerator(model.numberCutGenerators()-1)->setGlobalCuts(true); |
---|
1120 | } |
---|
1121 | } |
---|
1122 | // Do search |
---|
1123 | if (logLevel > 1) |
---|
1124 | model_->messageHandler()->message(CBC_START_SUB, model_->messages()) |
---|
1125 | << name |
---|
1126 | << model.getMaximumNodes() |
---|
1127 | << CoinMessageEol; |
---|
1128 | // probably faster to use a basis to get integer solutions |
---|
1129 | model.setSpecialOptions(model.specialOptions() | 2); |
---|
1130 | #ifdef CBC_THREAD |
---|
1131 | if (model_->getNumberThreads() > 0 && (model_->getThreadMode()&4) != 0) { |
---|
1132 | // See if at root node |
---|
1133 | bool atRoot = model_->getNodeCount() == 0; |
---|
1134 | int passNumber = model_->getCurrentPassNumber(); |
---|
1135 | if (atRoot && passNumber == 1) |
---|
1136 | model.setNumberThreads(model_->getNumberThreads()); |
---|
1137 | } |
---|
1138 | #endif |
---|
1139 | model.setParentModel(*model_); |
---|
1140 | model.setMaximumSolutions(maximumSolutions); |
---|
1141 | model.setOriginalColumns(process.originalColumns()); |
---|
1142 | model.setSearchStrategy(-1); |
---|
1143 | // If no feasibility pump then insert a lightweight one |
---|
1144 | if (feasibilityPumpOptions_ >= 0 || feasibilityPumpOptions_ == -2) { |
---|
1145 | CbcHeuristicFPump * fpump = NULL; |
---|
1146 | for (int i = 0; i < model.numberHeuristics(); i++) { |
---|
1147 | CbcHeuristicFPump* pump = |
---|
1148 | dynamic_cast<CbcHeuristicFPump*>(model.heuristic(i)); |
---|
1149 | if (pump) { |
---|
1150 | fpump = pump; |
---|
1151 | break; |
---|
1152 | } |
---|
1153 | } |
---|
1154 | if (!fpump) { |
---|
1155 | CbcHeuristicFPump heuristic4; |
---|
1156 | // use any cutoff |
---|
1157 | heuristic4.setFakeCutoff(0.5*COIN_DBL_MAX); |
---|
1158 | if (fractionSmall_<=1.0) |
---|
1159 | heuristic4.setMaximumPasses(10); |
---|
1160 | int pumpTune = feasibilityPumpOptions_; |
---|
1161 | if (pumpTune==-2) |
---|
1162 | pumpTune = 4; // proximity |
---|
1163 | if (pumpTune > 0) { |
---|
1164 | /* |
---|
1165 | >=10000000 for using obj |
---|
1166 | >=1000000 use as accumulate switch |
---|
1167 | >=1000 use index+1 as number of large loops |
---|
1168 | >=100 use 0.05 objvalue as increment |
---|
1169 | %100 == 10,20 etc for experimentation |
---|
1170 | 1 == fix ints at bounds, 2 fix all integral ints, 3 and continuous at bounds |
---|
1171 | 4 and static continuous, 5 as 3 but no internal integers |
---|
1172 | 6 as 3 but all slack basis! |
---|
1173 | */ |
---|
1174 | double value = solver2->getObjSense() * solver2->getObjValue(); |
---|
1175 | int w = pumpTune / 10; |
---|
1176 | int ix = w % 10; |
---|
1177 | w /= 10; |
---|
1178 | int c = w % 10; |
---|
1179 | w /= 10; |
---|
1180 | int r = w; |
---|
1181 | int accumulate = r / 1000; |
---|
1182 | r -= 1000 * accumulate; |
---|
1183 | if (accumulate >= 10) { |
---|
1184 | int which = accumulate / 10; |
---|
1185 | accumulate -= 10 * which; |
---|
1186 | which--; |
---|
1187 | // weights and factors |
---|
1188 | double weight[] = {0.1, 0.1, 0.5, 0.5, 1.0, 1.0, 5.0, 5.0}; |
---|
1189 | double factor[] = {0.1, 0.5, 0.1, 0.5, 0.1, 0.5, 0.1, 0.5}; |
---|
1190 | heuristic4.setInitialWeight(weight[which]); |
---|
1191 | heuristic4.setWeightFactor(factor[which]); |
---|
1192 | } |
---|
1193 | // fake cutoff |
---|
1194 | if (c) { |
---|
1195 | double cutoff; |
---|
1196 | solver2->getDblParam(OsiDualObjectiveLimit, cutoff); |
---|
1197 | cutoff = CoinMin(cutoff, value + 0.1 * fabs(value) * c); |
---|
1198 | heuristic4.setFakeCutoff(cutoff); |
---|
1199 | } |
---|
1200 | if (r) { |
---|
1201 | // also set increment |
---|
1202 | //double increment = (0.01*i+0.005)*(fabs(value)+1.0e-12); |
---|
1203 | double increment = 0.0; |
---|
1204 | heuristic4.setAbsoluteIncrement(increment); |
---|
1205 | heuristic4.setAccumulate(accumulate); |
---|
1206 | heuristic4.setMaximumRetries(r + 1); |
---|
1207 | } |
---|
1208 | pumpTune = pumpTune % 100; |
---|
1209 | if (pumpTune == 6) |
---|
1210 | pumpTune = 13; |
---|
1211 | if (pumpTune != 13) |
---|
1212 | pumpTune = pumpTune % 10; |
---|
1213 | heuristic4.setWhen(pumpTune); |
---|
1214 | if (ix) { |
---|
1215 | heuristic4.setFeasibilityPumpOptions(ix*10); |
---|
1216 | } |
---|
1217 | } |
---|
1218 | model.addHeuristic(&heuristic4, "feasibility pump", 0); |
---|
1219 | |
---|
1220 | } |
---|
1221 | } else if (feasibilityPumpOptions_==-3) { |
---|
1222 | // add all (except this) |
---|
1223 | for (int i = 0; i < model_->numberHeuristics(); i++) { |
---|
1224 | if (strcmp(heuristicName(),model_->heuristic(i)->heuristicName())) |
---|
1225 | model.addHeuristic(model_->heuristic(i)); |
---|
1226 | } |
---|
1227 | } |
---|
1228 | // modify heuristics |
---|
1229 | for (int i = 0; i < model.numberHeuristics(); i++) { |
---|
1230 | // reset lastNode |
---|
1231 | CbcHeuristicRINS * rins = |
---|
1232 | dynamic_cast<CbcHeuristicRINS*>(model.heuristic(i)); |
---|
1233 | if (rins) { |
---|
1234 | rins->setLastNode(-1000); |
---|
1235 | rins->setSolutionCount(0); |
---|
1236 | } |
---|
1237 | } |
---|
1238 | //printf("sol %x\n",inputSolution_); |
---|
1239 | #ifdef CGL_DEBUG |
---|
1240 | if ((model_->specialOptions()&1) != 0) { |
---|
1241 | const OsiRowCutDebugger *debugger = model.solver()->getRowCutDebugger() ; |
---|
1242 | if (debugger) { |
---|
1243 | printf("On optimal path CC\n"); |
---|
1244 | } |
---|
1245 | } |
---|
1246 | #endif |
---|
1247 | if (inputSolution_) { |
---|
1248 | // translate and add a serendipity heuristic |
---|
1249 | int numberColumns = solver2->getNumCols(); |
---|
1250 | const int * which = process.originalColumns(); |
---|
1251 | OsiSolverInterface * solver3 = solver2->clone(); |
---|
1252 | for (int i = 0; i < numberColumns; i++) { |
---|
1253 | if (isHeuristicInteger(solver3,i)) { |
---|
1254 | int k = which[i]; |
---|
1255 | double value = inputSolution_[k]; |
---|
1256 | //if (value) |
---|
1257 | //printf("orig col %d now %d val %g\n", |
---|
1258 | // k,i,value); |
---|
1259 | solver3->setColLower(i, value); |
---|
1260 | solver3->setColUpper(i, value); |
---|
1261 | } |
---|
1262 | } |
---|
1263 | solver3->setDblParam(OsiDualObjectiveLimit, COIN_DBL_MAX); |
---|
1264 | solver3->resolve(); |
---|
1265 | if (!solver3->isProvenOptimal()) { |
---|
1266 | // Try just setting nonzeros |
---|
1267 | OsiSolverInterface * solver4 = solver2->clone(); |
---|
1268 | for (int i = 0; i < numberColumns; i++) { |
---|
1269 | if (isHeuristicInteger(solver4,i)) { |
---|
1270 | int k = which[i]; |
---|
1271 | double value = floor(inputSolution_[k] + 0.5); |
---|
1272 | if (value) { |
---|
1273 | solver3->setColLower(i, value); |
---|
1274 | solver3->setColUpper(i, value); |
---|
1275 | } |
---|
1276 | } |
---|
1277 | } |
---|
1278 | solver4->setDblParam(OsiDualObjectiveLimit, COIN_DBL_MAX); |
---|
1279 | solver4->resolve(); |
---|
1280 | int nBad = -1; |
---|
1281 | if (solver4->isProvenOptimal()) { |
---|
1282 | nBad = 0; |
---|
1283 | const double * solution = solver4->getColSolution(); |
---|
1284 | for (int i = 0; i < numberColumns; i++) { |
---|
1285 | if (isHeuristicInteger(solver4,i)) { |
---|
1286 | double value = floor(solution[i] + 0.5); |
---|
1287 | if (fabs(value - solution[i]) > 1.0e-6) |
---|
1288 | nBad++; |
---|
1289 | } |
---|
1290 | } |
---|
1291 | } |
---|
1292 | if (nBad) { |
---|
1293 | delete solver4; |
---|
1294 | } else { |
---|
1295 | delete solver3; |
---|
1296 | solver3 = solver4; |
---|
1297 | } |
---|
1298 | } |
---|
1299 | if (solver3->isProvenOptimal()) { |
---|
1300 | // good |
---|
1301 | CbcSerendipity heuristic(model); |
---|
1302 | double value = solver3->getObjSense() * solver3->getObjValue(); |
---|
1303 | heuristic.setInputSolution(solver3->getColSolution(), value); |
---|
1304 | value = value + 1.0e-7*(1.0 + fabs(value)); |
---|
1305 | value *= solver3->getObjSense(); |
---|
1306 | model.setCutoff(value); |
---|
1307 | model.addHeuristic(&heuristic, "Previous solution", 0); |
---|
1308 | //printf("added seren\n"); |
---|
1309 | } else { |
---|
1310 | double value = model_->getMinimizationObjValue(); |
---|
1311 | value = value + 1.0e-7*(1.0 + fabs(value)); |
---|
1312 | value *= solver3->getObjSense(); |
---|
1313 | model.setCutoff(value); |
---|
1314 | sprintf(generalPrint, "Unable to insert previous solution - using cutoff of %g", |
---|
1315 | value); |
---|
1316 | model_->messageHandler()->message(CBC_FPUMP1, model_->messages()) |
---|
1317 | << generalPrint |
---|
1318 | << CoinMessageEol; |
---|
1319 | #ifdef CLP_INVESTIGATE |
---|
1320 | printf("NOT added seren\n"); |
---|
1321 | solver3->writeMps("bad_seren"); |
---|
1322 | solver->writeMps("orig_seren"); |
---|
1323 | #endif |
---|
1324 | } |
---|
1325 | delete solver3; |
---|
1326 | } |
---|
1327 | if (model_->searchStrategy() == 2) { |
---|
1328 | model.setNumberStrong(5); |
---|
1329 | model.setNumberBeforeTrust(5); |
---|
1330 | } |
---|
1331 | if (model.getNumCols()) { |
---|
1332 | if (numberNodes >= 0) { |
---|
1333 | setCutAndHeuristicOptions(model); |
---|
1334 | // not too many iterations |
---|
1335 | model.setMaximumNumberIterations(iterationMultiplier*(numberNodes + 10)); |
---|
1336 | // Not fast stuff |
---|
1337 | model.setFastNodeDepth(-1); |
---|
1338 | //model.solver()->writeMps("before"); |
---|
1339 | } else if (model.fastNodeDepth() >= 1000000) { |
---|
1340 | // already set |
---|
1341 | model.setFastNodeDepth(model.fastNodeDepth() - 1000000); |
---|
1342 | } |
---|
1343 | model.setWhenCuts(999998); |
---|
1344 | #define ALWAYS_DUAL |
---|
1345 | #ifdef ALWAYS_DUAL |
---|
1346 | OsiSolverInterface * solverD = model.solver(); |
---|
1347 | bool takeHint; |
---|
1348 | OsiHintStrength strength; |
---|
1349 | solverD->getHintParam(OsiDoDualInResolve, takeHint, strength); |
---|
1350 | solverD->setHintParam(OsiDoDualInResolve, true, OsiHintDo); |
---|
1351 | #endif |
---|
1352 | model.passInEventHandler(model_->getEventHandler()); |
---|
1353 | // say model_ is sitting there |
---|
1354 | int saveOptions = model_->specialOptions(); |
---|
1355 | model_->setSpecialOptions(saveOptions|1048576); |
---|
1356 | // and switch off debugger |
---|
1357 | model.setSpecialOptions(model.specialOptions()&(~1)); |
---|
1358 | #if 0 //def COIN_HAS_CLP |
---|
1359 | OsiClpSolverInterface * clpSolver |
---|
1360 | = dynamic_cast<OsiClpSolverInterface *> (model.solver()); |
---|
1361 | if (clpSolver) |
---|
1362 | clpSolver->zapDebugger(); |
---|
1363 | #endif |
---|
1364 | #ifdef CONFLICT_CUTS |
---|
1365 | if ((model_->moreSpecialOptions()&4194304)!=0) |
---|
1366 | model.zapGlobalCuts(); |
---|
1367 | #endif |
---|
1368 | #ifdef CGL_DEBUG |
---|
1369 | if ((model_->specialOptions()&1) != 0) { |
---|
1370 | const OsiRowCutDebugger *debugger = model.solver()->getRowCutDebugger() ; |
---|
1371 | if (debugger) { |
---|
1372 | printf("On optimal path DD\n"); |
---|
1373 | } |
---|
1374 | } |
---|
1375 | #endif |
---|
1376 | model.branchAndBound(); |
---|
1377 | model_->setHeuristicModel(NULL); |
---|
1378 | model_->setSpecialOptions(saveOptions); |
---|
1379 | #ifdef ALWAYS_DUAL |
---|
1380 | solverD = model.solver(); |
---|
1381 | solverD->setHintParam(OsiDoDualInResolve, takeHint, strength); |
---|
1382 | #endif |
---|
1383 | numberNodesDone_ = model.getNodeCount(); |
---|
1384 | #ifdef COIN_DEVELOP |
---|
1385 | printf("sub branch %d nodes, %d iterations - max %d\n", |
---|
1386 | model.getNodeCount(), model.getIterationCount(), |
---|
1387 | 100*(numberNodes + 10)); |
---|
1388 | #endif |
---|
1389 | if (numberNodes < 0) { |
---|
1390 | model_->incrementIterationCount(model.getIterationCount()); |
---|
1391 | model_->incrementNodeCount(model.getNodeCount()); |
---|
1392 | // update best solution (in case ctrl-c) |
---|
1393 | // !!! not a good idea - think a bit harder |
---|
1394 | //model_->setMinimizationObjValue(model.getMinimizationObjValue()); |
---|
1395 | for (int iGenerator = 0; iGenerator < model.numberCutGenerators(); iGenerator++) { |
---|
1396 | CbcCutGenerator * generator = model.cutGenerator(iGenerator); |
---|
1397 | sprintf(generalPrint, |
---|
1398 | "%s was tried %d times and created %d cuts of which %d were active after adding rounds of cuts (%.3f seconds)", |
---|
1399 | generator->cutGeneratorName(), |
---|
1400 | generator->numberTimesEntered(), |
---|
1401 | generator->numberCutsInTotal() + |
---|
1402 | generator->numberColumnCuts(), |
---|
1403 | generator->numberCutsActive(), |
---|
1404 | generator->timeInCutGenerator()); |
---|
1405 | CglStored * stored = dynamic_cast<CglStored*>(generator->generator()); |
---|
1406 | if (stored && !generator->numberCutsInTotal()) |
---|
1407 | continue; |
---|
1408 | #ifndef CLP_INVESTIGATE |
---|
1409 | CglImplication * implication = dynamic_cast<CglImplication*>(generator->generator()); |
---|
1410 | if (implication) |
---|
1411 | continue; |
---|
1412 | #endif |
---|
1413 | model_->messageHandler()->message(CBC_FPUMP1, model_->messages()) |
---|
1414 | << generalPrint |
---|
1415 | << CoinMessageEol; |
---|
1416 | } |
---|
1417 | } |
---|
1418 | } else { |
---|
1419 | // empty model |
---|
1420 | model.setMinimizationObjValue(model.solver()->getObjSense()*model.solver()->getObjValue()); |
---|
1421 | } |
---|
1422 | if (logLevel > 1) |
---|
1423 | model_->messageHandler()->message(CBC_END_SUB, model_->messages()) |
---|
1424 | << name |
---|
1425 | << CoinMessageEol; |
---|
1426 | if (model.getMinimizationObjValue() < CoinMin(cutoff, 1.0e30)) { |
---|
1427 | // solution |
---|
1428 | if (model.getNumCols()) |
---|
1429 | returnCode = model.isProvenOptimal() ? 3 : 1; |
---|
1430 | else |
---|
1431 | returnCode = 3; |
---|
1432 | // post process |
---|
1433 | #ifdef COIN_HAS_CLP |
---|
1434 | OsiClpSolverInterface * clpSolver = dynamic_cast< OsiClpSolverInterface*> (model.solver()); |
---|
1435 | if (clpSolver) { |
---|
1436 | ClpSimplex * lpSolver = clpSolver->getModelPtr(); |
---|
1437 | lpSolver->setSpecialOptions(lpSolver->specialOptions() | 0x01000000); // say is Cbc (and in branch and bound) |
---|
1438 | } |
---|
1439 | #endif |
---|
1440 | //if (fractionSmall_ < 1000000.0) |
---|
1441 | process.postProcess(*model.solver()); |
---|
1442 | if (solver->isProvenOptimal() && solver->getObjValue()*solver->getObjSense() < cutoff) { |
---|
1443 | // Solution now back in solver |
---|
1444 | int numberColumns = solver->getNumCols(); |
---|
1445 | memcpy(newSolution, solver->getColSolution(), |
---|
1446 | numberColumns*sizeof(double)); |
---|
1447 | newSolutionValue = model.getMinimizationObjValue(); |
---|
1448 | } else { |
---|
1449 | // odd - but no good |
---|
1450 | returnCode = 0; // so will be infeasible |
---|
1451 | } |
---|
1452 | } else { |
---|
1453 | // no good |
---|
1454 | returnCode = model.isProvenInfeasible() ? 2 : 0; // so will be infeasible |
---|
1455 | } |
---|
1456 | int totalNumberIterations = model.getIterationCount() + |
---|
1457 | process.numberIterationsPre() + |
---|
1458 | process.numberIterationsPost(); |
---|
1459 | if (totalNumberIterations > 100*(numberNodes + 10) |
---|
1460 | && fractionSmall_ < 1000000.0) { |
---|
1461 | // only allow smaller problems |
---|
1462 | fractionSmall = fractionSmall_; |
---|
1463 | fractionSmall_ *= 0.9; |
---|
1464 | #ifdef CLP_INVESTIGATE |
---|
1465 | printf("changing fractionSmall from %g to %g for %s as %d iterations\n", |
---|
1466 | fractionSmall, fractionSmall_, name.c_str(), totalNumberIterations); |
---|
1467 | #endif |
---|
1468 | } |
---|
1469 | if (model.status() == 5) |
---|
1470 | model_->sayEventHappened(); |
---|
1471 | #ifdef COIN_DEVELOP |
---|
1472 | if (model.isProvenInfeasible()) |
---|
1473 | status = 1; |
---|
1474 | else if (model.isProvenOptimal()) |
---|
1475 | status = 2; |
---|
1476 | #endif |
---|
1477 | } |
---|
1478 | } |
---|
1479 | } else { |
---|
1480 | returnCode = 2; // infeasible finished |
---|
1481 | if (logLevel > 1){ |
---|
1482 | printf("Infeasible on initial solve\n"); |
---|
1483 | } |
---|
1484 | } |
---|
1485 | model_->setSpecialOptions(saveModelOptions); |
---|
1486 | model_->setLogLevel(logLevel); |
---|
1487 | if (returnCode == 1 || returnCode == 2) { |
---|
1488 | OsiSolverInterface * solverC = model_->continuousSolver(); |
---|
1489 | if (false && solverC) { |
---|
1490 | const double * lower = solver->getColLower(); |
---|
1491 | const double * upper = solver->getColUpper(); |
---|
1492 | const double * lowerC = solverC->getColLower(); |
---|
1493 | const double * upperC = solverC->getColUpper(); |
---|
1494 | bool good = true; |
---|
1495 | for (int iColumn = 0; iColumn < numberColumns; iColumn++) { |
---|
1496 | if (isHeuristicInteger(solverC,iColumn)) { |
---|
1497 | if (lower[iColumn] > lowerC[iColumn] && |
---|
1498 | upper[iColumn] < upperC[iColumn]) { |
---|
1499 | good = false; |
---|
1500 | printf("CUT - can't add\n"); |
---|
1501 | break; |
---|
1502 | } |
---|
1503 | } |
---|
1504 | } |
---|
1505 | if (good) { |
---|
1506 | double * cut = new double [numberColumns]; |
---|
1507 | int * which = new int [numberColumns]; |
---|
1508 | double rhs = -1.0; |
---|
1509 | int n = 0; |
---|
1510 | for (int iColumn = 0; iColumn < numberColumns; iColumn++) { |
---|
1511 | if (isHeuristicInteger(solverC,iColumn)) { |
---|
1512 | if (lower[iColumn] == upperC[iColumn]) { |
---|
1513 | rhs += lower[iColumn]; |
---|
1514 | cut[n] = 1.0; |
---|
1515 | which[n++] = iColumn; |
---|
1516 | } else if (upper[iColumn] == lowerC[iColumn]) { |
---|
1517 | rhs -= upper[iColumn]; |
---|
1518 | cut[n] = -1.0; |
---|
1519 | which[n++] = iColumn; |
---|
1520 | } |
---|
1521 | } |
---|
1522 | } |
---|
1523 | printf("CUT has %d entries\n", n); |
---|
1524 | OsiRowCut newCut; |
---|
1525 | newCut.setLb(-COIN_DBL_MAX); |
---|
1526 | newCut.setUb(rhs); |
---|
1527 | newCut.setRow(n, which, cut, false); |
---|
1528 | model_->makeGlobalCut(newCut); |
---|
1529 | delete [] cut; |
---|
1530 | delete [] which; |
---|
1531 | } |
---|
1532 | } |
---|
1533 | #ifdef COIN_DEVELOP |
---|
1534 | if (status == 1) |
---|
1535 | printf("heuristic could add cut because infeasible (%s)\n", heuristicName_.c_str()); |
---|
1536 | else if (status == 2) |
---|
1537 | printf("heuristic could add cut because optimal (%s)\n", heuristicName_.c_str()); |
---|
1538 | #endif |
---|
1539 | } |
---|
1540 | if (reset) { |
---|
1541 | for (int iColumn = 0; iColumn < numberColumns; iColumn++) { |
---|
1542 | if (reset[iColumn]) |
---|
1543 | solver->setColLower(iColumn, 0.0); |
---|
1544 | } |
---|
1545 | delete [] reset; |
---|
1546 | } |
---|
1547 | #ifdef HISTORY_STATISTICS |
---|
1548 | getHistoryStatistics_ = true; |
---|
1549 | #endif |
---|
1550 | solver->setHintParam(OsiDoReducePrint, takeHint, strength); |
---|
1551 | return returnCode; |
---|
1552 | } |
---|
1553 | // Set input solution |
---|
1554 | void |
---|
1555 | CbcHeuristic::setInputSolution(const double * solution, double objValue) |
---|
1556 | { |
---|
1557 | delete [] inputSolution_; |
---|
1558 | inputSolution_ = NULL; |
---|
1559 | if (model_ && solution) { |
---|
1560 | int numberColumns = model_->getNumCols(); |
---|
1561 | inputSolution_ = new double [numberColumns+1]; |
---|
1562 | memcpy(inputSolution_, solution, numberColumns*sizeof(double)); |
---|
1563 | inputSolution_[numberColumns] = objValue; |
---|
1564 | } |
---|
1565 | } |
---|
1566 | |
---|
1567 | //############################################################################## |
---|
1568 | |
---|
1569 | inline int compare3BranchingObjects(const CbcBranchingObject* br0, |
---|
1570 | const CbcBranchingObject* br1) |
---|
1571 | { |
---|
1572 | const int t0 = br0->type(); |
---|
1573 | const int t1 = br1->type(); |
---|
1574 | if (t0 < t1) { |
---|
1575 | return -1; |
---|
1576 | } |
---|
1577 | if (t0 > t1) { |
---|
1578 | return 1; |
---|
1579 | } |
---|
1580 | return br0->compareOriginalObject(br1); |
---|
1581 | } |
---|
1582 | |
---|
1583 | //============================================================================== |
---|
1584 | |
---|
1585 | inline bool compareBranchingObjects(const CbcBranchingObject* br0, |
---|
1586 | const CbcBranchingObject* br1) |
---|
1587 | { |
---|
1588 | return compare3BranchingObjects(br0, br1) < 0; |
---|
1589 | } |
---|
1590 | |
---|
1591 | //============================================================================== |
---|
1592 | |
---|
1593 | void |
---|
1594 | CbcHeuristicNode::gutsOfConstructor(CbcModel& model) |
---|
1595 | { |
---|
1596 | // CbcHeurDebugNodes(&model); |
---|
1597 | CbcNode* node = model.currentNode(); |
---|
1598 | brObj_ = new CbcBranchingObject*[node->depth()]; |
---|
1599 | CbcNodeInfo* nodeInfo = node->nodeInfo(); |
---|
1600 | int cnt = 0; |
---|
1601 | while (nodeInfo->parentBranch() != NULL) { |
---|
1602 | const OsiBranchingObject* br = nodeInfo->parentBranch(); |
---|
1603 | const CbcBranchingObject* cbcbr = dynamic_cast<const CbcBranchingObject*>(br); |
---|
1604 | if (! cbcbr) { |
---|
1605 | throw CoinError("CbcHeuristicNode can be used only with CbcBranchingObjects.\n", |
---|
1606 | "gutsOfConstructor", |
---|
1607 | "CbcHeuristicNode", |
---|
1608 | __FILE__, __LINE__); |
---|
1609 | } |
---|
1610 | brObj_[cnt] = cbcbr->clone(); |
---|
1611 | brObj_[cnt]->previousBranch(); |
---|
1612 | ++cnt; |
---|
1613 | nodeInfo = nodeInfo->parent(); |
---|
1614 | } |
---|
1615 | std::sort(brObj_, brObj_ + cnt, compareBranchingObjects); |
---|
1616 | if (cnt <= 1) { |
---|
1617 | numObjects_ = cnt; |
---|
1618 | } else { |
---|
1619 | numObjects_ = 0; |
---|
1620 | CbcBranchingObject* br = NULL; // What should this be? |
---|
1621 | for (int i = 1; i < cnt; ++i) { |
---|
1622 | if (compare3BranchingObjects(brObj_[numObjects_], brObj_[i]) == 0) { |
---|
1623 | int comp = brObj_[numObjects_]->compareBranchingObject(brObj_[i], br != 0); |
---|
1624 | switch (comp) { |
---|
1625 | case CbcRangeSame: // the same range |
---|
1626 | case CbcRangeDisjoint: // disjoint decisions |
---|
1627 | // should not happen! we are on a chain! |
---|
1628 | abort(); |
---|
1629 | case CbcRangeSubset: // brObj_[numObjects_] is a subset of brObj_[i] |
---|
1630 | delete brObj_[i]; |
---|
1631 | break; |
---|
1632 | case CbcRangeSuperset: // brObj_[i] is a subset of brObj_[numObjects_] |
---|
1633 | delete brObj_[numObjects_]; |
---|
1634 | brObj_[numObjects_] = brObj_[i]; |
---|
1635 | break; |
---|
1636 | case CbcRangeOverlap: // overlap |
---|
1637 | delete brObj_[i]; |
---|
1638 | delete brObj_[numObjects_]; |
---|
1639 | brObj_[numObjects_] = br; |
---|
1640 | break; |
---|
1641 | } |
---|
1642 | continue; |
---|
1643 | } else { |
---|
1644 | brObj_[++numObjects_] = brObj_[i]; |
---|
1645 | } |
---|
1646 | } |
---|
1647 | ++numObjects_; |
---|
1648 | } |
---|
1649 | } |
---|
1650 | |
---|
1651 | //============================================================================== |
---|
1652 | |
---|
1653 | CbcHeuristicNode::CbcHeuristicNode(CbcModel& model) |
---|
1654 | { |
---|
1655 | gutsOfConstructor(model); |
---|
1656 | } |
---|
1657 | |
---|
1658 | //============================================================================== |
---|
1659 | |
---|
1660 | double |
---|
1661 | CbcHeuristicNode::distance(const CbcHeuristicNode* node) const |
---|
1662 | { |
---|
1663 | |
---|
1664 | const double disjointWeight = 1; |
---|
1665 | const double overlapWeight = 0.4; |
---|
1666 | const double subsetWeight = 0.2; |
---|
1667 | int countDisjointWeight = 0; |
---|
1668 | int countOverlapWeight = 0; |
---|
1669 | int countSubsetWeight = 0; |
---|
1670 | int i = 0; |
---|
1671 | int j = 0; |
---|
1672 | double dist = 0.0; |
---|
1673 | #ifdef PRINT_DEBUG |
---|
1674 | printf(" numObjects_ = %i, node->numObjects_ = %i\n", |
---|
1675 | numObjects_, node->numObjects_); |
---|
1676 | #endif |
---|
1677 | while ( i < numObjects_ && j < node->numObjects_) { |
---|
1678 | CbcBranchingObject* br0 = brObj_[i]; |
---|
1679 | const CbcBranchingObject* br1 = node->brObj_[j]; |
---|
1680 | #ifdef PRINT_DEBUG |
---|
1681 | const CbcIntegerBranchingObject* brPrint0 = |
---|
1682 | dynamic_cast<const CbcIntegerBranchingObject*>(br0); |
---|
1683 | const double* downBounds = brPrint0->downBounds(); |
---|
1684 | const double* upBounds = brPrint0->upBounds(); |
---|
1685 | int variable = brPrint0->variable(); |
---|
1686 | int way = brPrint0->way(); |
---|
1687 | printf(" br0: var %i downBd [%i,%i] upBd [%i,%i] way %i\n", |
---|
1688 | variable, static_cast<int>(downBounds[0]), static_cast<int>(downBounds[1]), |
---|
1689 | static_cast<int>(upBounds[0]), static_cast<int>(upBounds[1]), way); |
---|
1690 | const CbcIntegerBranchingObject* brPrint1 = |
---|
1691 | dynamic_cast<const CbcIntegerBranchingObject*>(br1); |
---|
1692 | downBounds = brPrint1->downBounds(); |
---|
1693 | upBounds = brPrint1->upBounds(); |
---|
1694 | variable = brPrint1->variable(); |
---|
1695 | way = brPrint1->way(); |
---|
1696 | printf(" br1: var %i downBd [%i,%i] upBd [%i,%i] way %i\n", |
---|
1697 | variable, static_cast<int>(downBounds[0]), static_cast<int>(downBounds[1]), |
---|
1698 | static_cast<int>(upBounds[0]), static_cast<int>(upBounds[1]), way); |
---|
1699 | #endif |
---|
1700 | const int brComp = compare3BranchingObjects(br0, br1); |
---|
1701 | if (brComp < 0) { |
---|
1702 | dist += subsetWeight; |
---|
1703 | countSubsetWeight++; |
---|
1704 | ++i; |
---|
1705 | } else if (brComp > 0) { |
---|
1706 | dist += subsetWeight; |
---|
1707 | countSubsetWeight++; |
---|
1708 | ++j; |
---|
1709 | } else { |
---|
1710 | const int comp = br0->compareBranchingObject(br1, false); |
---|
1711 | switch (comp) { |
---|
1712 | case CbcRangeSame: |
---|
1713 | // do nothing |
---|
1714 | break; |
---|
1715 | case CbcRangeDisjoint: // disjoint decisions |
---|
1716 | dist += disjointWeight; |
---|
1717 | countDisjointWeight++; |
---|
1718 | break; |
---|
1719 | case CbcRangeSubset: // subset one way or another |
---|
1720 | case CbcRangeSuperset: |
---|
1721 | dist += subsetWeight; |
---|
1722 | countSubsetWeight++; |
---|
1723 | break; |
---|
1724 | case CbcRangeOverlap: // overlap |
---|
1725 | dist += overlapWeight; |
---|
1726 | countOverlapWeight++; |
---|
1727 | break; |
---|
1728 | } |
---|
1729 | ++i; |
---|
1730 | ++j; |
---|
1731 | } |
---|
1732 | } |
---|
1733 | dist += subsetWeight * (numObjects_ - i + node->numObjects_ - j); |
---|
1734 | countSubsetWeight += (numObjects_ - i + node->numObjects_ - j); |
---|
1735 | COIN_DETAIL_PRINT(printf("subset = %i, overlap = %i, disjoint = %i\n", countSubsetWeight, |
---|
1736 | countOverlapWeight, countDisjointWeight)); |
---|
1737 | return dist; |
---|
1738 | } |
---|
1739 | |
---|
1740 | //============================================================================== |
---|
1741 | |
---|
1742 | CbcHeuristicNode::~CbcHeuristicNode() |
---|
1743 | { |
---|
1744 | for (int i = 0; i < numObjects_; ++i) { |
---|
1745 | delete brObj_[i]; |
---|
1746 | } |
---|
1747 | delete [] brObj_; |
---|
1748 | } |
---|
1749 | |
---|
1750 | //============================================================================== |
---|
1751 | |
---|
1752 | double |
---|
1753 | CbcHeuristicNode::minDistance(const CbcHeuristicNodeList& nodeList) const |
---|
1754 | { |
---|
1755 | double minDist = COIN_DBL_MAX; |
---|
1756 | for (int i = nodeList.size() - 1; i >= 0; --i) { |
---|
1757 | minDist = CoinMin(minDist, distance(nodeList.node(i))); |
---|
1758 | } |
---|
1759 | return minDist; |
---|
1760 | } |
---|
1761 | |
---|
1762 | //============================================================================== |
---|
1763 | |
---|
1764 | bool |
---|
1765 | CbcHeuristicNode::minDistanceIsSmall(const CbcHeuristicNodeList& nodeList, |
---|
1766 | const double threshold) const |
---|
1767 | { |
---|
1768 | for (int i = nodeList.size() - 1; i >= 0; --i) { |
---|
1769 | if (distance(nodeList.node(i)) >= threshold) { |
---|
1770 | continue; |
---|
1771 | } else { |
---|
1772 | return true; |
---|
1773 | } |
---|
1774 | } |
---|
1775 | return false; |
---|
1776 | } |
---|
1777 | |
---|
1778 | //============================================================================== |
---|
1779 | |
---|
1780 | double |
---|
1781 | CbcHeuristicNode::avgDistance(const CbcHeuristicNodeList& nodeList) const |
---|
1782 | { |
---|
1783 | if (nodeList.size() == 0) { |
---|
1784 | return COIN_DBL_MAX; |
---|
1785 | } |
---|
1786 | double sumDist = 0; |
---|
1787 | for (int i = nodeList.size() - 1; i >= 0; --i) { |
---|
1788 | sumDist += distance(nodeList.node(i)); |
---|
1789 | } |
---|
1790 | return sumDist / nodeList.size(); |
---|
1791 | } |
---|
1792 | |
---|
1793 | //############################################################################## |
---|
1794 | |
---|
1795 | // Default Constructor |
---|
1796 | CbcRounding::CbcRounding() |
---|
1797 | : CbcHeuristic() |
---|
1798 | { |
---|
1799 | // matrix and row copy will automatically be empty |
---|
1800 | seed_ = 7654321; |
---|
1801 | down_ = NULL; |
---|
1802 | up_ = NULL; |
---|
1803 | equal_ = NULL; |
---|
1804 | //whereFrom_ |= 16*(1+256); // allow more often |
---|
1805 | } |
---|
1806 | |
---|
1807 | // Constructor from model |
---|
1808 | CbcRounding::CbcRounding(CbcModel & model) |
---|
1809 | : CbcHeuristic(model) |
---|
1810 | { |
---|
1811 | // Get a copy of original matrix (and by row for rounding); |
---|
1812 | assert(model.solver()); |
---|
1813 | if (model.solver()->getNumRows()) { |
---|
1814 | matrix_ = *model.solver()->getMatrixByCol(); |
---|
1815 | matrixByRow_ = *model.solver()->getMatrixByRow(); |
---|
1816 | validate(); |
---|
1817 | } |
---|
1818 | down_ = NULL; |
---|
1819 | up_ = NULL; |
---|
1820 | equal_ = NULL; |
---|
1821 | seed_ = 7654321; |
---|
1822 | //whereFrom_ |= 16*(1+256); // allow more often |
---|
1823 | } |
---|
1824 | |
---|
1825 | // Destructor |
---|
1826 | CbcRounding::~CbcRounding () |
---|
1827 | { |
---|
1828 | delete [] down_; |
---|
1829 | delete [] up_; |
---|
1830 | delete [] equal_; |
---|
1831 | } |
---|
1832 | |
---|
1833 | // Clone |
---|
1834 | CbcHeuristic * |
---|
1835 | CbcRounding::clone() const |
---|
1836 | { |
---|
1837 | return new CbcRounding(*this); |
---|
1838 | } |
---|
1839 | // Create C++ lines to get to current state |
---|
1840 | void |
---|
1841 | CbcRounding::generateCpp( FILE * fp) |
---|
1842 | { |
---|
1843 | CbcRounding other; |
---|
1844 | fprintf(fp, "0#include \"CbcHeuristic.hpp\"\n"); |
---|
1845 | fprintf(fp, "3 CbcRounding rounding(*cbcModel);\n"); |
---|
1846 | CbcHeuristic::generateCpp(fp, "rounding"); |
---|
1847 | if (seed_ != other.seed_) |
---|
1848 | fprintf(fp, "3 rounding.setSeed(%d);\n", seed_); |
---|
1849 | else |
---|
1850 | fprintf(fp, "4 rounding.setSeed(%d);\n", seed_); |
---|
1851 | fprintf(fp, "3 cbcModel->addHeuristic(&rounding);\n"); |
---|
1852 | } |
---|
1853 | //#define NEW_ROUNDING |
---|
1854 | // Copy constructor |
---|
1855 | CbcRounding::CbcRounding(const CbcRounding & rhs) |
---|
1856 | : |
---|
1857 | CbcHeuristic(rhs), |
---|
1858 | matrix_(rhs.matrix_), |
---|
1859 | matrixByRow_(rhs.matrixByRow_), |
---|
1860 | seed_(rhs.seed_) |
---|
1861 | { |
---|
1862 | #ifdef NEW_ROUNDING |
---|
1863 | int numberColumns = matrix_.getNumCols(); |
---|
1864 | down_ = CoinCopyOfArray(rhs.down_, numberColumns); |
---|
1865 | up_ = CoinCopyOfArray(rhs.up_, numberColumns); |
---|
1866 | equal_ = CoinCopyOfArray(rhs.equal_, numberColumns); |
---|
1867 | #else |
---|
1868 | down_ = NULL; |
---|
1869 | up_ = NULL; |
---|
1870 | equal_ = NULL; |
---|
1871 | #endif |
---|
1872 | } |
---|
1873 | |
---|
1874 | // Assignment operator |
---|
1875 | CbcRounding & |
---|
1876 | CbcRounding::operator=( const CbcRounding & rhs) |
---|
1877 | { |
---|
1878 | if (this != &rhs) { |
---|
1879 | CbcHeuristic::operator=(rhs); |
---|
1880 | matrix_ = rhs.matrix_; |
---|
1881 | matrixByRow_ = rhs.matrixByRow_; |
---|
1882 | #ifdef NEW_ROUNDING |
---|
1883 | delete [] down_; |
---|
1884 | delete [] up_; |
---|
1885 | delete [] equal_; |
---|
1886 | int numberColumns = matrix_.getNumCols(); |
---|
1887 | down_ = CoinCopyOfArray(rhs.down_, numberColumns); |
---|
1888 | up_ = CoinCopyOfArray(rhs.up_, numberColumns); |
---|
1889 | equal_ = CoinCopyOfArray(rhs.equal_, numberColumns); |
---|
1890 | #else |
---|
1891 | down_ = NULL; |
---|
1892 | up_ = NULL; |
---|
1893 | equal_ = NULL; |
---|
1894 | #endif |
---|
1895 | seed_ = rhs.seed_; |
---|
1896 | } |
---|
1897 | return *this; |
---|
1898 | } |
---|
1899 | |
---|
1900 | // Resets stuff if model changes |
---|
1901 | void |
---|
1902 | CbcRounding::resetModel(CbcModel * model) |
---|
1903 | { |
---|
1904 | model_ = model; |
---|
1905 | // Get a copy of original matrix (and by row for rounding); |
---|
1906 | assert(model_->solver()); |
---|
1907 | matrix_ = *model_->solver()->getMatrixByCol(); |
---|
1908 | matrixByRow_ = *model_->solver()->getMatrixByRow(); |
---|
1909 | validate(); |
---|
1910 | } |
---|
1911 | /* Check whether the heuristic should run at all |
---|
1912 | 0 - before cuts at root node (or from doHeuristics) |
---|
1913 | 1 - during cuts at root |
---|
1914 | 2 - after root node cuts |
---|
1915 | 3 - after cuts at other nodes |
---|
1916 | 4 - during cuts at other nodes |
---|
1917 | 8 added if previous heuristic in loop found solution |
---|
1918 | */ |
---|
1919 | bool |
---|
1920 | CbcRounding::shouldHeurRun(int whereFrom) |
---|
1921 | { |
---|
1922 | if (whereFrom!=4) { |
---|
1923 | return CbcHeuristic::shouldHeurRun(whereFrom); |
---|
1924 | } else { |
---|
1925 | numCouldRun_++; |
---|
1926 | return shouldHeurRun_randomChoice(); |
---|
1927 | } |
---|
1928 | } |
---|
1929 | // See if rounding will give solution |
---|
1930 | // Sets value of solution |
---|
1931 | // Assumes rhs for original matrix still okay |
---|
1932 | // At present only works with integers |
---|
1933 | // Fix values if asked for |
---|
1934 | // Returns 1 if solution, 0 if not |
---|
1935 | int |
---|
1936 | CbcRounding::solution(double & solutionValue, |
---|
1937 | double * betterSolution) |
---|
1938 | { |
---|
1939 | |
---|
1940 | numCouldRun_++; |
---|
1941 | // See if to do |
---|
1942 | if (!when() || (when() % 10 == 1 && model_->phase() != 1) || |
---|
1943 | (when() % 10 == 2 && (model_->phase() != 2 && model_->phase() != 3))) |
---|
1944 | return 0; // switched off |
---|
1945 | numRuns_++; |
---|
1946 | #ifdef HEURISTIC_INFORM |
---|
1947 | printf("Entering heuristic %s - nRuns %d numCould %d when %d\n", |
---|
1948 | heuristicName(),numRuns_,numCouldRun_,when_); |
---|
1949 | #endif |
---|
1950 | OsiSolverInterface * solver = model_->solver(); |
---|
1951 | double direction = solver->getObjSense(); |
---|
1952 | double newSolutionValue = direction * solver->getObjValue(); |
---|
1953 | return solution(solutionValue, betterSolution, newSolutionValue); |
---|
1954 | } |
---|
1955 | // See if rounding will give solution |
---|
1956 | // Sets value of solution |
---|
1957 | // Assumes rhs for original matrix still okay |
---|
1958 | // At present only works with integers |
---|
1959 | // Fix values if asked for |
---|
1960 | // Returns 1 if solution, 0 if not |
---|
1961 | int |
---|
1962 | CbcRounding::solution(double & solutionValue, |
---|
1963 | double * betterSolution, |
---|
1964 | double newSolutionValue) |
---|
1965 | { |
---|
1966 | |
---|
1967 | // See if to do |
---|
1968 | if (!when() || (when() % 10 == 1 && model_->phase() != 1) || |
---|
1969 | (when() % 10 == 2 && (model_->phase() != 2 && model_->phase() != 3))) |
---|
1970 | return 0; // switched off |
---|
1971 | OsiSolverInterface * solver = model_->solver(); |
---|
1972 | const double * lower = solver->getColLower(); |
---|
1973 | const double * upper = solver->getColUpper(); |
---|
1974 | const double * rowLower = solver->getRowLower(); |
---|
1975 | const double * rowUpper = solver->getRowUpper(); |
---|
1976 | const double * solution = solver->getColSolution(); |
---|
1977 | const double * objective = solver->getObjCoefficients(); |
---|
1978 | double integerTolerance = model_->getDblParam(CbcModel::CbcIntegerTolerance); |
---|
1979 | double primalTolerance; |
---|
1980 | solver->getDblParam(OsiPrimalTolerance, primalTolerance); |
---|
1981 | double useTolerance = primalTolerance; |
---|
1982 | |
---|
1983 | int numberRows = matrix_.getNumRows(); |
---|
1984 | assert (numberRows <= solver->getNumRows()); |
---|
1985 | if (numberRows == 0){ |
---|
1986 | return 0; |
---|
1987 | } |
---|
1988 | int numberIntegers = model_->numberIntegers(); |
---|
1989 | const int * integerVariable = model_->integerVariable(); |
---|
1990 | int i; |
---|
1991 | double direction = solver->getObjSense(); |
---|
1992 | //double newSolutionValue = direction*solver->getObjValue(); |
---|
1993 | int returnCode = 0; |
---|
1994 | // Column copy |
---|
1995 | const double * element = matrix_.getElements(); |
---|
1996 | const int * row = matrix_.getIndices(); |
---|
1997 | const CoinBigIndex * columnStart = matrix_.getVectorStarts(); |
---|
1998 | const int * columnLength = matrix_.getVectorLengths(); |
---|
1999 | // Row copy |
---|
2000 | const double * elementByRow = matrixByRow_.getElements(); |
---|
2001 | const int * column = matrixByRow_.getIndices(); |
---|
2002 | const CoinBigIndex * rowStart = matrixByRow_.getVectorStarts(); |
---|
2003 | const int * rowLength = matrixByRow_.getVectorLengths(); |
---|
2004 | |
---|
2005 | // Get solution array for heuristic solution |
---|
2006 | int numberColumns = solver->getNumCols(); |
---|
2007 | double * newSolution = new double [numberColumns]; |
---|
2008 | memcpy(newSolution, solution, numberColumns*sizeof(double)); |
---|
2009 | |
---|
2010 | double * rowActivity = new double[numberRows]; |
---|
2011 | memset(rowActivity, 0, numberRows*sizeof(double)); |
---|
2012 | for (i = 0; i < numberColumns; i++) { |
---|
2013 | int j; |
---|
2014 | double value = newSolution[i]; |
---|
2015 | if (value < lower[i]) { |
---|
2016 | value = lower[i]; |
---|
2017 | newSolution[i] = value; |
---|
2018 | } else if (value > upper[i]) { |
---|
2019 | value = upper[i]; |
---|
2020 | newSolution[i] = value; |
---|
2021 | } |
---|
2022 | if (value) { |
---|
2023 | for (j = columnStart[i]; |
---|
2024 | j < columnStart[i] + columnLength[i]; j++) { |
---|
2025 | int iRow = row[j]; |
---|
2026 | rowActivity[iRow] += value * element[j]; |
---|
2027 | } |
---|
2028 | } |
---|
2029 | } |
---|
2030 | // check was feasible - if not adjust (cleaning may move) |
---|
2031 | for (i = 0; i < numberRows; i++) { |
---|
2032 | if (rowActivity[i] < rowLower[i]) { |
---|
2033 | //assert (rowActivity[i]>rowLower[i]-1000.0*primalTolerance); |
---|
2034 | rowActivity[i] = rowLower[i]; |
---|
2035 | } else if (rowActivity[i] > rowUpper[i]) { |
---|
2036 | //assert (rowActivity[i]<rowUpper[i]+1000.0*primalTolerance); |
---|
2037 | rowActivity[i] = rowUpper[i]; |
---|
2038 | } |
---|
2039 | } |
---|
2040 | for (i = 0; i < numberIntegers; i++) { |
---|
2041 | int iColumn = integerVariable[i]; |
---|
2042 | double value = newSolution[iColumn]; |
---|
2043 | //double thisTolerance = integerTolerance; |
---|
2044 | if (fabs(floor(value + 0.5) - value) > integerTolerance) { |
---|
2045 | double below = floor(value); |
---|
2046 | double newValue = newSolution[iColumn]; |
---|
2047 | double cost = direction * objective[iColumn]; |
---|
2048 | double move; |
---|
2049 | if (cost > 0.0) { |
---|
2050 | // try up |
---|
2051 | move = 1.0 - (value - below); |
---|
2052 | } else if (cost < 0.0) { |
---|
2053 | // try down |
---|
2054 | move = below - value; |
---|
2055 | } else { |
---|
2056 | // won't be able to move unless we can grab another variable |
---|
2057 | double randomNumber = randomNumberGenerator_.randomDouble(); |
---|
2058 | // which way? |
---|
2059 | if (randomNumber < 0.5) |
---|
2060 | move = below - value; |
---|
2061 | else |
---|
2062 | move = 1.0 - (value - below); |
---|
2063 | } |
---|
2064 | newValue += move; |
---|
2065 | newSolution[iColumn] = newValue; |
---|
2066 | newSolutionValue += move * cost; |
---|
2067 | int j; |
---|
2068 | for (j = columnStart[iColumn]; |
---|
2069 | j < columnStart[iColumn] + columnLength[iColumn]; j++) { |
---|
2070 | int iRow = row[j]; |
---|
2071 | rowActivity[iRow] += move * element[j]; |
---|
2072 | } |
---|
2073 | } |
---|
2074 | } |
---|
2075 | |
---|
2076 | double penalty = 0.0; |
---|
2077 | // see if feasible - just using singletons |
---|
2078 | for (i = 0; i < numberRows; i++) { |
---|
2079 | double value = rowActivity[i]; |
---|
2080 | double thisInfeasibility = 0.0; |
---|
2081 | if (value < rowLower[i] - primalTolerance) |
---|
2082 | thisInfeasibility = value - rowLower[i]; |
---|
2083 | else if (value > rowUpper[i] + primalTolerance) |
---|
2084 | thisInfeasibility = value - rowUpper[i]; |
---|
2085 | if (thisInfeasibility) { |
---|
2086 | // See if there are any slacks I can use to fix up |
---|
2087 | // maybe put in coding for multiple slacks? |
---|
2088 | double bestCost = 1.0e50; |
---|
2089 | int k; |
---|
2090 | int iBest = -1; |
---|
2091 | double addCost = 0.0; |
---|
2092 | double newValue = 0.0; |
---|
2093 | double changeRowActivity = 0.0; |
---|
2094 | double absInfeasibility = fabs(thisInfeasibility); |
---|
2095 | for (k = rowStart[i]; k < rowStart[i] + rowLength[i]; k++) { |
---|
2096 | int iColumn = column[k]; |
---|
2097 | // See if all elements help |
---|
2098 | if (columnLength[iColumn] == 1) { |
---|
2099 | double currentValue = newSolution[iColumn]; |
---|
2100 | double elementValue = elementByRow[k]; |
---|
2101 | double lowerValue = lower[iColumn]; |
---|
2102 | double upperValue = upper[iColumn]; |
---|
2103 | double gap = rowUpper[i] - rowLower[i]; |
---|
2104 | double absElement = fabs(elementValue); |
---|
2105 | if (thisInfeasibility*elementValue > 0.0) { |
---|
2106 | // we want to reduce |
---|
2107 | if ((currentValue - lowerValue)*absElement >= absInfeasibility) { |
---|
2108 | // possible - check if integer |
---|
2109 | double distance = absInfeasibility / absElement; |
---|
2110 | double thisCost = -direction * objective[iColumn] * distance; |
---|
2111 | if (isHeuristicInteger(solver,iColumn)) { |
---|
2112 | distance = ceil(distance - useTolerance); |
---|
2113 | if (currentValue - distance >= lowerValue - useTolerance) { |
---|
2114 | if (absInfeasibility - distance*absElement < -gap - useTolerance) |
---|
2115 | thisCost = 1.0e100; // no good |
---|
2116 | else |
---|
2117 | thisCost = -direction * objective[iColumn] * distance; |
---|
2118 | } else { |
---|
2119 | thisCost = 1.0e100; // no good |
---|
2120 | } |
---|
2121 | } |
---|
2122 | if (thisCost < bestCost) { |
---|
2123 | bestCost = thisCost; |
---|
2124 | iBest = iColumn; |
---|
2125 | addCost = thisCost; |
---|
2126 | newValue = currentValue - distance; |
---|
2127 | changeRowActivity = -distance * elementValue; |
---|
2128 | } |
---|
2129 | } |
---|
2130 | } else { |
---|
2131 | // we want to increase |
---|
2132 | if ((upperValue - currentValue)*absElement >= absInfeasibility) { |
---|
2133 | // possible - check if integer |
---|
2134 | double distance = absInfeasibility / absElement; |
---|
2135 | double thisCost = direction * objective[iColumn] * distance; |
---|
2136 | if (isHeuristicInteger(solver,iColumn)) { |
---|
2137 | distance = ceil(distance - 1.0e-7); |
---|
2138 | assert (currentValue - distance <= upperValue + useTolerance); |
---|
2139 | if (absInfeasibility - distance*absElement < -gap - useTolerance) |
---|
2140 | thisCost = 1.0e100; // no good |
---|
2141 | else |
---|
2142 | thisCost = direction * objective[iColumn] * distance; |
---|
2143 | } |
---|
2144 | if (thisCost < bestCost) { |
---|
2145 | bestCost = thisCost; |
---|
2146 | iBest = iColumn; |
---|
2147 | addCost = thisCost; |
---|
2148 | newValue = currentValue + distance; |
---|
2149 | changeRowActivity = distance * elementValue; |
---|
2150 | } |
---|
2151 | } |
---|
2152 | } |
---|
2153 | } |
---|
2154 | } |
---|
2155 | if (iBest >= 0) { |
---|
2156 | /*printf("Infeasibility of %g on row %d cost %g\n", |
---|
2157 | thisInfeasibility,i,addCost);*/ |
---|
2158 | newSolution[iBest] = newValue; |
---|
2159 | thisInfeasibility = 0.0; |
---|
2160 | newSolutionValue += addCost; |
---|
2161 | rowActivity[i] += changeRowActivity; |
---|
2162 | } |
---|
2163 | penalty += fabs(thisInfeasibility); |
---|
2164 | } |
---|
2165 | } |
---|
2166 | if (penalty) { |
---|
2167 | // see if feasible using any |
---|
2168 | // first continuous |
---|
2169 | double penaltyChange = 0.0; |
---|
2170 | int iColumn; |
---|
2171 | for (iColumn = 0; iColumn < numberColumns; iColumn++) { |
---|
2172 | if (isHeuristicInteger(solver,iColumn)) |
---|
2173 | continue; |
---|
2174 | double currentValue = newSolution[iColumn]; |
---|
2175 | double lowerValue = lower[iColumn]; |
---|
2176 | double upperValue = upper[iColumn]; |
---|
2177 | int j; |
---|
2178 | int anyBadDown = 0; |
---|
2179 | int anyBadUp = 0; |
---|
2180 | double upImprovement = 0.0; |
---|
2181 | double downImprovement = 0.0; |
---|
2182 | for (j = columnStart[iColumn]; |
---|
2183 | j < columnStart[iColumn] + columnLength[iColumn]; j++) { |
---|
2184 | int iRow = row[j]; |
---|
2185 | if (rowUpper[iRow] > rowLower[iRow]) { |
---|
2186 | double value = element[j]; |
---|
2187 | if (rowActivity[iRow] > rowUpper[iRow] + primalTolerance) { |
---|
2188 | // infeasible above |
---|
2189 | downImprovement += value; |
---|
2190 | upImprovement -= value; |
---|
2191 | if (value > 0.0) |
---|
2192 | anyBadUp++; |
---|
2193 | else |
---|
2194 | anyBadDown++; |
---|
2195 | } else if (rowActivity[iRow] > rowUpper[iRow] - primalTolerance) { |
---|
2196 | // feasible at ub |
---|
2197 | if (value > 0.0) { |
---|
2198 | upImprovement -= value; |
---|
2199 | anyBadUp++; |
---|
2200 | } else { |
---|
2201 | downImprovement += value; |
---|
2202 | anyBadDown++; |
---|
2203 | } |
---|
2204 | } else if (rowActivity[iRow] > rowLower[iRow] + primalTolerance) { |
---|
2205 | // feasible in interior |
---|
2206 | } else if (rowActivity[iRow] > rowLower[iRow] - primalTolerance) { |
---|
2207 | // feasible at lb |
---|
2208 | if (value < 0.0) { |
---|
2209 | upImprovement += value; |
---|
2210 | anyBadUp++; |
---|
2211 | } else { |
---|
2212 | downImprovement -= value; |
---|
2213 | anyBadDown++; |
---|
2214 | } |
---|
2215 | } else { |
---|
2216 | // infeasible below |
---|
2217 | downImprovement -= value; |
---|
2218 | upImprovement += value; |
---|
2219 | if (value < 0.0) |
---|
2220 | anyBadUp++; |
---|
2221 | else |
---|
2222 | anyBadDown++; |
---|
2223 | } |
---|
2224 | } else { |
---|
2225 | // equality row |
---|
2226 | double value = element[j]; |
---|
2227 | if (rowActivity[iRow] > rowUpper[iRow] + primalTolerance) { |
---|
2228 | // infeasible above |
---|
2229 | downImprovement += value; |
---|
2230 | upImprovement -= value; |
---|
2231 | if (value > 0.0) |
---|
2232 | anyBadUp++; |
---|
2233 | else |
---|
2234 | anyBadDown++; |
---|
2235 | } else if (rowActivity[iRow] < rowLower[iRow] - primalTolerance) { |
---|
2236 | // infeasible below |
---|
2237 | downImprovement -= value; |
---|
2238 | upImprovement += value; |
---|
2239 | if (value < 0.0) |
---|
2240 | anyBadUp++; |
---|
2241 | else |
---|
2242 | anyBadDown++; |
---|
2243 | } else { |
---|
2244 | // feasible - no good |
---|
2245 | anyBadUp = -1; |
---|
2246 | anyBadDown = -1; |
---|
2247 | break; |
---|
2248 | } |
---|
2249 | } |
---|
2250 | } |
---|
2251 | // could change tests for anyBad |
---|
2252 | if (anyBadUp) |
---|
2253 | upImprovement = 0.0; |
---|
2254 | if (anyBadDown) |
---|
2255 | downImprovement = 0.0; |
---|
2256 | double way = 0.0; |
---|
2257 | double improvement = 0.0; |
---|
2258 | if (downImprovement > 0.0 && currentValue > lowerValue) { |
---|
2259 | way = -1.0; |
---|
2260 | improvement = downImprovement; |
---|
2261 | } else if (upImprovement > 0.0 && currentValue < upperValue) { |
---|
2262 | way = 1.0; |
---|
2263 | improvement = upImprovement; |
---|
2264 | } |
---|
2265 | if (way) { |
---|
2266 | // can improve |
---|
2267 | double distance; |
---|
2268 | if (way > 0.0) |
---|
2269 | distance = upperValue - currentValue; |
---|
2270 | else |
---|
2271 | distance = currentValue - lowerValue; |
---|
2272 | for (j = columnStart[iColumn]; |
---|
2273 | j < columnStart[iColumn] + columnLength[iColumn]; j++) { |
---|
2274 | int iRow = row[j]; |
---|
2275 | double value = element[j] * way; |
---|
2276 | if (rowActivity[iRow] > rowUpper[iRow] + primalTolerance) { |
---|
2277 | // infeasible above |
---|
2278 | assert (value < 0.0); |
---|
2279 | double gap = rowActivity[iRow] - rowUpper[iRow]; |
---|
2280 | if (gap + value*distance < 0.0) |
---|
2281 | distance = -gap / value; |
---|
2282 | } else if (rowActivity[iRow] < rowLower[iRow] - primalTolerance) { |
---|
2283 | // infeasible below |
---|
2284 | assert (value > 0.0); |
---|
2285 | double gap = rowActivity[iRow] - rowLower[iRow]; |
---|
2286 | if (gap + value*distance > 0.0) |
---|
2287 | distance = -gap / value; |
---|
2288 | } else { |
---|
2289 | // feasible |
---|
2290 | if (value > 0) { |
---|
2291 | double gap = rowActivity[iRow] - rowUpper[iRow]; |
---|
2292 | if (gap + value*distance > 0.0) |
---|
2293 | distance = -gap / value; |
---|
2294 | } else { |
---|
2295 | double gap = rowActivity[iRow] - rowLower[iRow]; |
---|
2296 | if (gap + value*distance < 0.0) |
---|
2297 | distance = -gap / value; |
---|
2298 | } |
---|
2299 | } |
---|
2300 | } |
---|
2301 | //move |
---|
2302 | penaltyChange += improvement * distance; |
---|
2303 | distance *= way; |
---|
2304 | newSolution[iColumn] += distance; |
---|
2305 | newSolutionValue += direction * objective[iColumn] * distance; |
---|
2306 | for (j = columnStart[iColumn]; |
---|
2307 | j < columnStart[iColumn] + columnLength[iColumn]; j++) { |
---|
2308 | int iRow = row[j]; |
---|
2309 | double value = element[j]; |
---|
2310 | rowActivity[iRow] += distance * value; |
---|
2311 | } |
---|
2312 | } |
---|
2313 | } |
---|
2314 | // and now all if improving |
---|
2315 | double lastChange = penaltyChange ? 1.0 : 0.0; |
---|
2316 | int numberPasses=0; |
---|
2317 | while (lastChange > 1.0e-2 && numberPasses < 1000) { |
---|
2318 | lastChange = 0; |
---|
2319 | numberPasses++; |
---|
2320 | for (iColumn = 0; iColumn < numberColumns; iColumn++) { |
---|
2321 | bool isInteger = isHeuristicInteger(solver,iColumn); |
---|
2322 | double currentValue = newSolution[iColumn]; |
---|
2323 | double lowerValue = lower[iColumn]; |
---|
2324 | double upperValue = upper[iColumn]; |
---|
2325 | int j; |
---|
2326 | int anyBadDown = 0; |
---|
2327 | int anyBadUp = 0; |
---|
2328 | double upImprovement = 0.0; |
---|
2329 | double downImprovement = 0.0; |
---|
2330 | for (j = columnStart[iColumn]; |
---|
2331 | j < columnStart[iColumn] + columnLength[iColumn]; j++) { |
---|
2332 | int iRow = row[j]; |
---|
2333 | double value = element[j]; |
---|
2334 | if (isInteger) { |
---|
2335 | if (value > 0.0) { |
---|
2336 | if (rowActivity[iRow] + value > rowUpper[iRow] + primalTolerance) |
---|
2337 | anyBadUp++; |
---|
2338 | if (rowActivity[iRow] - value < rowLower[iRow] - primalTolerance) |
---|
2339 | anyBadDown++; |
---|
2340 | } else { |
---|
2341 | if (rowActivity[iRow] - value > rowUpper[iRow] + primalTolerance) |
---|
2342 | anyBadDown++; |
---|
2343 | if (rowActivity[iRow] + value < rowLower[iRow] - primalTolerance) |
---|
2344 | anyBadUp++; |
---|
2345 | } |
---|
2346 | } |
---|
2347 | if (rowUpper[iRow] > rowLower[iRow]) { |
---|
2348 | if (rowActivity[iRow] > rowUpper[iRow] + primalTolerance) { |
---|
2349 | // infeasible above |
---|
2350 | downImprovement += value; |
---|
2351 | upImprovement -= value; |
---|
2352 | if (value > 0.0) |
---|
2353 | anyBadUp++; |
---|
2354 | else |
---|
2355 | anyBadDown++; |
---|
2356 | } else if (rowActivity[iRow] > rowUpper[iRow] - primalTolerance) { |
---|
2357 | // feasible at ub |
---|
2358 | if (value > 0.0) { |
---|
2359 | upImprovement -= value; |
---|
2360 | anyBadUp++; |
---|
2361 | } else { |
---|
2362 | downImprovement += value; |
---|
2363 | anyBadDown++; |
---|
2364 | } |
---|
2365 | } else if (rowActivity[iRow] > rowLower[iRow] + primalTolerance) { |
---|
2366 | // feasible in interior |
---|
2367 | } else if (rowActivity[iRow] > rowLower[iRow] - primalTolerance) { |
---|
2368 | // feasible at lb |
---|
2369 | if (value < 0.0) { |
---|
2370 | upImprovement += value; |
---|
2371 | anyBadUp++; |
---|
2372 | } else { |
---|
2373 | downImprovement -= value; |
---|
2374 | anyBadDown++; |
---|
2375 | } |
---|
2376 | } else { |
---|
2377 | // infeasible below |
---|
2378 | downImprovement -= value; |
---|
2379 | upImprovement += value; |
---|
2380 | if (value < 0.0) |
---|
2381 | anyBadUp++; |
---|
2382 | else |
---|
2383 | anyBadDown++; |
---|
2384 | } |
---|
2385 | } else { |
---|
2386 | // equality row |
---|
2387 | if (rowActivity[iRow] > rowUpper[iRow] + primalTolerance) { |
---|
2388 | // infeasible above |
---|
2389 | downImprovement += value; |
---|
2390 | upImprovement -= value; |
---|
2391 | if (value > 0.0) |
---|
2392 | anyBadUp++; |
---|
2393 | else |
---|
2394 | anyBadDown++; |
---|
2395 | } else if (rowActivity[iRow] < rowLower[iRow] - primalTolerance) { |
---|
2396 | // infeasible below |
---|
2397 | downImprovement -= value; |
---|
2398 | upImprovement += value; |
---|
2399 | if (value < 0.0) |
---|
2400 | anyBadUp++; |
---|
2401 | else |
---|
2402 | anyBadDown++; |
---|
2403 | } else { |
---|
2404 | // feasible - no good |
---|
2405 | anyBadUp = -1; |
---|
2406 | anyBadDown = -1; |
---|
2407 | break; |
---|
2408 | } |
---|
2409 | } |
---|
2410 | } |
---|
2411 | // could change tests for anyBad |
---|
2412 | if (anyBadUp) |
---|
2413 | upImprovement = 0.0; |
---|
2414 | if (anyBadDown) |
---|
2415 | downImprovement = 0.0; |
---|
2416 | double way = 0.0; |
---|
2417 | double improvement = 0.0; |
---|
2418 | if (downImprovement > 0.0 && currentValue > lowerValue) { |
---|
2419 | way = -1.0; |
---|
2420 | improvement = downImprovement; |
---|
2421 | if (isInteger&¤tValue<lowerValue+0.99) |
---|
2422 | continue; // no good |
---|
2423 | } else if (upImprovement > 0.0 && currentValue < upperValue) { |
---|
2424 | way = 1.0; |
---|
2425 | improvement = upImprovement; |
---|
2426 | if (isInteger&¤tValue>upperValue-0.99) |
---|
2427 | continue; // no good |
---|
2428 | } |
---|
2429 | if (way) { |
---|
2430 | // can improve |
---|
2431 | double distance = COIN_DBL_MAX; |
---|
2432 | for (j = columnStart[iColumn]; |
---|
2433 | j < columnStart[iColumn] + columnLength[iColumn]; j++) { |
---|
2434 | int iRow = row[j]; |
---|
2435 | double value = element[j] * way; |
---|
2436 | if (rowActivity[iRow] > rowUpper[iRow] + primalTolerance) { |
---|
2437 | // infeasible above |
---|
2438 | assert (value < 0.0); |
---|
2439 | double gap = rowActivity[iRow] - rowUpper[iRow]; |
---|
2440 | if (gap + value*distance < 0.0) { |
---|
2441 | // If integer then has to move by 1 |
---|
2442 | if (!isInteger) |
---|
2443 | distance = -gap / value; |
---|
2444 | else |
---|
2445 | distance = CoinMax(-gap / value, 1.0); |
---|
2446 | } |
---|
2447 | } else if (rowActivity[iRow] < rowLower[iRow] - primalTolerance) { |
---|
2448 | // infeasible below |
---|
2449 | assert (value > 0.0); |
---|
2450 | double gap = rowActivity[iRow] - rowLower[iRow]; |
---|
2451 | if (gap + value*distance > 0.0) { |
---|
2452 | // If integer then has to move by 1 |
---|
2453 | if (!isInteger) |
---|
2454 | distance = -gap / value; |
---|
2455 | else |
---|
2456 | distance = CoinMax(-gap / value, 1.0); |
---|
2457 | } |
---|
2458 | } else { |
---|
2459 | // feasible |
---|
2460 | if (value > 0) { |
---|
2461 | double gap = rowActivity[iRow] - rowUpper[iRow]; |
---|
2462 | if (gap + value*distance > 0.0) |
---|
2463 | distance = -gap / value; |
---|
2464 | } else { |
---|
2465 | double gap = rowActivity[iRow] - rowLower[iRow]; |
---|
2466 | if (gap + value*distance < 0.0) |
---|
2467 | distance = -gap / value; |
---|
2468 | } |
---|
2469 | } |
---|
2470 | } |
---|
2471 | if (isInteger) |
---|
2472 | distance = floor(distance + 1.05e-8); |
---|
2473 | if (!distance) { |
---|
2474 | // should never happen |
---|
2475 | //printf("zero distance in CbcRounding - debug\n"); |
---|
2476 | } |
---|
2477 | //move |
---|
2478 | lastChange += improvement * distance; |
---|
2479 | distance *= way; |
---|
2480 | newSolution[iColumn] += distance; |
---|
2481 | newSolutionValue += direction * objective[iColumn] * distance; |
---|
2482 | for (j = columnStart[iColumn]; |
---|
2483 | j < columnStart[iColumn] + columnLength[iColumn]; j++) { |
---|
2484 | int iRow = row[j]; |
---|
2485 | double value = element[j]; |
---|
2486 | rowActivity[iRow] += distance * value; |
---|
2487 | } |
---|
2488 | } |
---|
2489 | } |
---|
2490 | penaltyChange += lastChange; |
---|
2491 | } |
---|
2492 | penalty -= penaltyChange; |
---|
2493 | if (penalty < 1.0e-5*fabs(penaltyChange)) { |
---|
2494 | // recompute |
---|
2495 | penalty = 0.0; |
---|
2496 | for (i = 0; i < numberRows; i++) { |
---|
2497 | double value = rowActivity[i]; |
---|
2498 | if (value < rowLower[i] - primalTolerance) |
---|
2499 | penalty += rowLower[i] - value; |
---|
2500 | else if (value > rowUpper[i] + primalTolerance) |
---|
2501 | penalty += value - rowUpper[i]; |
---|
2502 | } |
---|
2503 | } |
---|
2504 | } |
---|
2505 | |
---|
2506 | // Could also set SOS (using random) and repeat |
---|
2507 | if (!penalty) { |
---|
2508 | // See if we can do better |
---|
2509 | //seed_++; |
---|
2510 | //CoinSeedRandom(seed_); |
---|
2511 | // Random number between 0 and 1. |
---|
2512 | double randomNumber = randomNumberGenerator_.randomDouble(); |
---|
2513 | int iPass; |
---|
2514 | int start[2]; |
---|
2515 | int end[2]; |
---|
2516 | int iRandom = static_cast<int> (randomNumber * (static_cast<double> (numberIntegers))); |
---|
2517 | start[0] = iRandom; |
---|
2518 | end[0] = numberIntegers; |
---|
2519 | start[1] = 0; |
---|
2520 | end[1] = iRandom; |
---|
2521 | for (iPass = 0; iPass < 2; iPass++) { |
---|
2522 | int i; |
---|
2523 | for (i = start[iPass]; i < end[iPass]; i++) { |
---|
2524 | int iColumn = integerVariable[i]; |
---|
2525 | #ifndef NDEBUG |
---|
2526 | double value = newSolution[iColumn]; |
---|
2527 | assert (fabs(floor(value + 0.5) - value) < integerTolerance); |
---|
2528 | #endif |
---|
2529 | double cost = direction * objective[iColumn]; |
---|
2530 | double move = 0.0; |
---|
2531 | if (cost > 0.0) |
---|
2532 | move = -1.0; |
---|
2533 | else if (cost < 0.0) |
---|
2534 | move = 1.0; |
---|
2535 | while (move) { |
---|
2536 | bool good = true; |
---|
2537 | double newValue = newSolution[iColumn] + move; |
---|
2538 | if (newValue < lower[iColumn] - useTolerance || |
---|
2539 | newValue > upper[iColumn] + useTolerance) { |
---|
2540 | move = 0.0; |
---|
2541 | } else { |
---|
2542 | // see if we can move |
---|
2543 | int j; |
---|
2544 | for (j = columnStart[iColumn]; |
---|
2545 | j < columnStart[iColumn] + columnLength[iColumn]; j++) { |
---|
2546 | int iRow = row[j]; |
---|
2547 | double newActivity = rowActivity[iRow] + move * element[j]; |
---|
2548 | if (newActivity < rowLower[iRow] - primalTolerance || |
---|
2549 | newActivity > rowUpper[iRow] + primalTolerance) { |
---|
2550 | good = false; |
---|
2551 | break; |
---|
2552 | } |
---|
2553 | } |
---|
2554 | if (good) { |
---|
2555 | newSolution[iColumn] = newValue; |
---|
2556 | newSolutionValue += move * cost; |
---|
2557 | int j; |
---|
2558 | for (j = columnStart[iColumn]; |
---|
2559 | j < columnStart[iColumn] + columnLength[iColumn]; j++) { |
---|
2560 | int iRow = row[j]; |
---|
2561 | rowActivity[iRow] += move * element[j]; |
---|
2562 | } |
---|
2563 | } else { |
---|
2564 | move = 0.0; |
---|
2565 | } |
---|
2566 | } |
---|
2567 | } |
---|
2568 | } |
---|
2569 | } |
---|
2570 | // Just in case of some stupidity |
---|
2571 | double objOffset = 0.0; |
---|
2572 | solver->getDblParam(OsiObjOffset, objOffset); |
---|
2573 | newSolutionValue = -objOffset; |
---|
2574 | for ( i = 0 ; i < numberColumns ; i++ ) |
---|
2575 | newSolutionValue += objective[i] * newSolution[i]; |
---|
2576 | newSolutionValue *= direction; |
---|
2577 | //printf("new solution value %g %g\n",newSolutionValue,solutionValue); |
---|
2578 | if (newSolutionValue < solutionValue) { |
---|
2579 | // paranoid check |
---|
2580 | memset(rowActivity, 0, numberRows*sizeof(double)); |
---|
2581 | for (i = 0; i < numberColumns; i++) { |
---|
2582 | int j; |
---|
2583 | double value = newSolution[i]; |
---|
2584 | if (value) { |
---|
2585 | for (j = columnStart[i]; |
---|
2586 | j < columnStart[i] + columnLength[i]; j++) { |
---|
2587 | int iRow = row[j]; |
---|
2588 | rowActivity[iRow] += value * element[j]; |
---|
2589 | } |
---|
2590 | } |
---|
2591 | } |
---|
2592 | // check was approximately feasible |
---|
2593 | bool feasible = true; |
---|
2594 | for (i = 0; i < numberRows; i++) { |
---|
2595 | if (rowActivity[i] < rowLower[i]) { |
---|
2596 | if (rowActivity[i] < rowLower[i] - 1000.0*primalTolerance) |
---|
2597 | feasible = false; |
---|
2598 | } else if (rowActivity[i] > rowUpper[i]) { |
---|
2599 | if (rowActivity[i] > rowUpper[i] + 1000.0*primalTolerance) |
---|
2600 | feasible = false; |
---|
2601 | } |
---|
2602 | } |
---|
2603 | if (feasible) { |
---|
2604 | // new solution |
---|
2605 | memcpy(betterSolution, newSolution, numberColumns*sizeof(double)); |
---|
2606 | solutionValue = newSolutionValue; |
---|
2607 | //printf("** Solution of %g found by rounding\n",newSolutionValue); |
---|
2608 | returnCode = 1; |
---|
2609 | } else { |
---|
2610 | // Can easily happen |
---|
2611 | //printf("Debug CbcRounding giving bad solution\n"); |
---|
2612 | } |
---|
2613 | } |
---|
2614 | } |
---|
2615 | #ifdef NEW_ROUNDING |
---|
2616 | if (!returnCode) { |
---|
2617 | #ifdef JJF_ZERO |
---|
2618 | // back to starting point |
---|
2619 | memcpy(newSolution, solution, numberColumns*sizeof(double)); |
---|
2620 | memset(rowActivity, 0, numberRows*sizeof(double)); |
---|
2621 | for (i = 0; i < numberColumns; i++) { |
---|
2622 | int j; |
---|
2623 | double value = newSolution[i]; |
---|
2624 | if (value < lower[i]) { |
---|
2625 | value = lower[i]; |
---|
2626 | newSolution[i] = value; |
---|
2627 | } else if (value > upper[i]) { |
---|
2628 | value = upper[i]; |
---|
2629 | newSolution[i] = value; |
---|
2630 | } |
---|
2631 | if (value) { |
---|
2632 | for (j = columnStart[i]; |
---|
2633 | j < columnStart[i] + columnLength[i]; j++) { |
---|
2634 | int iRow = row[j]; |
---|
2635 | rowActivity[iRow] += value * element[j]; |
---|
2636 | } |
---|
2637 | } |
---|
2638 | } |
---|
2639 | // check was feasible - if not adjust (cleaning may move) |
---|
2640 | for (i = 0; i < numberRows; i++) { |
---|
2641 | if (rowActivity[i] < rowLower[i]) { |
---|
2642 | //assert (rowActivity[i]>rowLower[i]-1000.0*primalTolerance); |
---|
2643 | rowActivity[i] = rowLower[i]; |
---|
2644 | } else if (rowActivity[i] > rowUpper[i]) { |
---|
2645 | //assert (rowActivity[i]<rowUpper[i]+1000.0*primalTolerance); |
---|
2646 | rowActivity[i] = rowUpper[i]; |
---|
2647 | } |
---|
2648 | } |
---|
2649 | #endif |
---|
2650 | int * candidate = new int [numberColumns]; |
---|
2651 | int nCandidate = 0; |
---|
2652 | for (iColumn = 0; iColumn < numberColumns; iColumn++) { |
---|
2653 | bool isInteger = isHeuristicInteger(solver,iColumn); |
---|
2654 | if (isInteger) { |
---|
2655 | double currentValue = newSolution[iColumn]; |
---|
2656 | if (fabs(currentValue - floor(currentValue + 0.5)) > 1.0e-8) |
---|
2657 | candidate[nCandidate++] = iColumn; |
---|
2658 | } |
---|
2659 | } |
---|
2660 | if (true) { |
---|
2661 | // Rounding as in Berthold |
---|
2662 | while (nCandidate) { |
---|
2663 | double infeasibility = 1.0e-7; |
---|
2664 | int iRow = -1; |
---|
2665 | for (i = 0; i < numberRows; i++) { |
---|
2666 | double value = 0.0; |
---|
2667 | if (rowActivity[i] < rowLower[i]) { |
---|
2668 | value = rowLower[i] - rowActivity[i]; |
---|
2669 | } else if (rowActivity[i] > rowUpper[i]) { |
---|
2670 | value = rowActivity[i] - rowUpper[i]; |
---|
2671 | } |
---|
2672 | if (value > infeasibility) { |
---|
2673 | infeasibility = value; |
---|
2674 | iRow = i; |
---|
2675 | } |
---|
2676 | } |
---|
2677 | if (iRow >= 0) { |
---|
2678 | // infeasible |
---|
2679 | } else { |
---|
2680 | // feasible |
---|
2681 | } |
---|
2682 | } |
---|
2683 | } else { |
---|
2684 | // Shifting as in Berthold |
---|
2685 | } |
---|
2686 | delete [] candidate; |
---|
2687 | } |
---|
2688 | #endif |
---|
2689 | delete [] newSolution; |
---|
2690 | delete [] rowActivity; |
---|
2691 | return returnCode; |
---|
2692 | } |
---|
2693 | // update model |
---|
2694 | void CbcRounding::setModel(CbcModel * model) |
---|
2695 | { |
---|
2696 | model_ = model; |
---|
2697 | // Get a copy of original matrix (and by row for rounding); |
---|
2698 | assert(model_->solver()); |
---|
2699 | if (model_->solver()->getNumRows()) { |
---|
2700 | matrix_ = *model_->solver()->getMatrixByCol(); |
---|
2701 | matrixByRow_ = *model_->solver()->getMatrixByRow(); |
---|
2702 | // make sure model okay for heuristic |
---|
2703 | validate(); |
---|
2704 | } |
---|
2705 | } |
---|
2706 | // Validate model i.e. sets when_ to 0 if necessary (may be NULL) |
---|
2707 | void |
---|
2708 | CbcRounding::validate() |
---|
2709 | { |
---|
2710 | if (model_ && (when() % 100) < 10) { |
---|
2711 | if (model_->numberIntegers() != |
---|
2712 | model_->numberObjects() && (model_->numberObjects() || |
---|
2713 | (model_->specialOptions()&1024) == 0)) { |
---|
2714 | int numberOdd = 0; |
---|
2715 | for (int i = 0; i < model_->numberObjects(); i++) { |
---|
2716 | if (!model_->object(i)->canDoHeuristics()) |
---|
2717 | numberOdd++; |
---|
2718 | } |
---|
2719 | if (numberOdd) |
---|
2720 | setWhen(0); |
---|
2721 | } |
---|
2722 | } |
---|
2723 | #ifdef NEW_ROUNDING |
---|
2724 | int numberColumns = matrix_.getNumCols(); |
---|
2725 | down_ = new unsigned short [numberColumns]; |
---|
2726 | up_ = new unsigned short [numberColumns]; |
---|
2727 | equal_ = new unsigned short [numberColumns]; |
---|
2728 | // Column copy |
---|
2729 | const double * element = matrix_.getElements(); |
---|
2730 | const int * row = matrix_.getIndices(); |
---|
2731 | const CoinBigIndex * columnStart = matrix_.getVectorStarts(); |
---|
2732 | const int * columnLength = matrix_.getVectorLengths(); |
---|
2733 | const double * rowLower = model.solver()->getRowLower(); |
---|
2734 | const double * rowUpper = model.solver()->getRowUpper(); |
---|
2735 | for (int i = 0; i < numberColumns; i++) { |
---|
2736 | int down = 0; |
---|
2737 | int up = 0; |
---|
2738 | int equal = 0; |
---|
2739 | if (columnLength[i] > 65535) { |
---|
2740 | equal[0] = 65535; |
---|
2741 | break; // unlikely to work |
---|
2742 | } |
---|
2743 | for (CoinBigIndex j = columnStart[i]; |
---|
2744 | j < columnStart[i] + columnLength[i]; j++) { |
---|
2745 | int iRow = row[j]; |
---|
2746 | if (rowLower[iRow] > -1.0e20 && rowUpper[iRow] < 1.0e20) { |
---|
2747 | equal++; |
---|
2748 | } else if (element[j] > 0.0) { |
---|
2749 | if (rowUpper[iRow] < 1.0e20) |
---|
2750 | up++; |
---|
2751 | else |
---|
2752 | down--; |
---|
2753 | } else { |
---|
2754 | if (rowLower[iRow] > -1.0e20) |
---|
2755 | up++; |
---|
2756 | else |
---|
2757 | down--; |
---|
2758 | } |
---|
2759 | } |
---|
2760 | down_[i] = (unsigned short) down; |
---|
2761 | up_[i] = (unsigned short) up; |
---|
2762 | equal_[i] = (unsigned short) equal; |
---|
2763 | } |
---|
2764 | #else |
---|
2765 | down_ = NULL; |
---|
2766 | up_ = NULL; |
---|
2767 | equal_ = NULL; |
---|
2768 | #endif |
---|
2769 | } |
---|
2770 | |
---|
2771 | // Default Constructor |
---|
2772 | CbcHeuristicPartial::CbcHeuristicPartial() |
---|
2773 | : CbcHeuristic() |
---|
2774 | { |
---|
2775 | fixPriority_ = 10000; |
---|
2776 | } |
---|
2777 | |
---|
2778 | // Constructor from model |
---|
2779 | CbcHeuristicPartial::CbcHeuristicPartial(CbcModel & model, int fixPriority, int numberNodes) |
---|
2780 | : CbcHeuristic(model) |
---|
2781 | { |
---|
2782 | fixPriority_ = fixPriority; |
---|
2783 | setNumberNodes(numberNodes); |
---|
2784 | validate(); |
---|
2785 | } |
---|
2786 | |
---|
2787 | // Destructor |
---|
2788 | CbcHeuristicPartial::~CbcHeuristicPartial () |
---|
2789 | { |
---|
2790 | } |
---|
2791 | |
---|
2792 | // Clone |
---|
2793 | CbcHeuristic * |
---|
2794 | CbcHeuristicPartial::clone() const |
---|
2795 | { |
---|
2796 | return new CbcHeuristicPartial(*this); |
---|
2797 | } |
---|
2798 | // Create C++ lines to get to current state |
---|
2799 | void |
---|
2800 | CbcHeuristicPartial::generateCpp( FILE * fp) |
---|
2801 | { |
---|
2802 | CbcHeuristicPartial other; |
---|
2803 | fprintf(fp, "0#include \"CbcHeuristic.hpp\"\n"); |
---|
2804 | fprintf(fp, "3 CbcHeuristicPartial partial(*cbcModel);\n"); |
---|
2805 | CbcHeuristic::generateCpp(fp, "partial"); |
---|
2806 | if (fixPriority_ != other.fixPriority_) |
---|
2807 | fprintf(fp, "3 partial.setFixPriority(%d);\n", fixPriority_); |
---|
2808 | else |
---|
2809 | fprintf(fp, "4 partial.setFixPriority(%d);\n", fixPriority_); |
---|
2810 | fprintf(fp, "3 cbcModel->addHeuristic(&partial);\n"); |
---|
2811 | } |
---|
2812 | //#define NEW_PARTIAL |
---|
2813 | // Copy constructor |
---|
2814 | CbcHeuristicPartial::CbcHeuristicPartial(const CbcHeuristicPartial & rhs) |
---|
2815 | : |
---|
2816 | CbcHeuristic(rhs), |
---|
2817 | fixPriority_(rhs.fixPriority_) |
---|
2818 | { |
---|
2819 | } |
---|
2820 | |
---|
2821 | // Assignment operator |
---|
2822 | CbcHeuristicPartial & |
---|
2823 | CbcHeuristicPartial::operator=( const CbcHeuristicPartial & rhs) |
---|
2824 | { |
---|
2825 | if (this != &rhs) { |
---|
2826 | CbcHeuristic::operator=(rhs); |
---|
2827 | fixPriority_ = rhs.fixPriority_; |
---|
2828 | } |
---|
2829 | return *this; |
---|
2830 | } |
---|
2831 | |
---|
2832 | // Resets stuff if model changes |
---|
2833 | void |
---|
2834 | CbcHeuristicPartial::resetModel(CbcModel * model) |
---|
2835 | { |
---|
2836 | model_ = model; |
---|
2837 | // Get a copy of original matrix (and by row for partial); |
---|
2838 | assert(model_->solver()); |
---|
2839 | validate(); |
---|
2840 | } |
---|
2841 | // See if partial will give solution |
---|
2842 | // Sets value of solution |
---|
2843 | // Assumes rhs for original matrix still okay |
---|
2844 | // At present only works with integers |
---|
2845 | // Fix values if asked for |
---|
2846 | // Returns 1 if solution, 0 if not |
---|
2847 | int |
---|
2848 | CbcHeuristicPartial::solution(double & solutionValue, |
---|
2849 | double * betterSolution) |
---|
2850 | { |
---|
2851 | // Return if already done |
---|
2852 | if (fixPriority_ < 0) |
---|
2853 | return 0; // switched off |
---|
2854 | #ifdef HEURISTIC_INFORM |
---|
2855 | printf("Entering heuristic %s - nRuns %d numCould %d when %d\n", |
---|
2856 | heuristicName(),numRuns_,numCouldRun_,when_); |
---|
2857 | #endif |
---|
2858 | const double * hotstartSolution = model_->hotstartSolution(); |
---|
2859 | const int * hotstartPriorities = model_->hotstartPriorities(); |
---|
2860 | if (!hotstartSolution) |
---|
2861 | return 0; |
---|
2862 | OsiSolverInterface * solver = model_->solver(); |
---|
2863 | |
---|
2864 | int numberIntegers = model_->numberIntegers(); |
---|
2865 | const int * integerVariable = model_->integerVariable(); |
---|
2866 | |
---|
2867 | OsiSolverInterface * newSolver = model_->continuousSolver()->clone(); |
---|
2868 | const double * colLower = newSolver->getColLower(); |
---|
2869 | const double * colUpper = newSolver->getColUpper(); |
---|
2870 | |
---|
2871 | double primalTolerance; |
---|
2872 | solver->getDblParam(OsiPrimalTolerance, primalTolerance); |
---|
2873 | |
---|
2874 | int i; |
---|
2875 | int numberFixed = 0; |
---|
2876 | int returnCode = 0; |
---|
2877 | |
---|
2878 | for (i = 0; i < numberIntegers; i++) { |
---|
2879 | int iColumn = integerVariable[i]; |
---|
2880 | if (abs(hotstartPriorities[iColumn]) <= fixPriority_) { |
---|
2881 | double value = hotstartSolution[iColumn]; |
---|
2882 | double lower = colLower[iColumn]; |
---|
2883 | double upper = colUpper[iColumn]; |
---|
2884 | value = CoinMax(value, lower); |
---|
2885 | value = CoinMin(value, upper); |
---|
2886 | if (fabs(value - floor(value + 0.5)) < 1.0e-8) { |
---|
2887 | value = floor(value + 0.5); |
---|
2888 | newSolver->setColLower(iColumn, value); |
---|
2889 | newSolver->setColUpper(iColumn, value); |
---|
2890 | numberFixed++; |
---|
2891 | } |
---|
2892 | } |
---|
2893 | } |
---|
2894 | if (numberFixed > numberIntegers / 5 - 100000000) { |
---|
2895 | #ifdef COIN_DEVELOP |
---|
2896 | printf("%d integers fixed\n", numberFixed); |
---|
2897 | #endif |
---|
2898 | returnCode = smallBranchAndBound(newSolver, numberNodes_, betterSolution, solutionValue, |
---|
2899 | model_->getCutoff(), "CbcHeuristicPartial"); |
---|
2900 | if (returnCode < 0) |
---|
2901 | returnCode = 0; // returned on size |
---|
2902 | //printf("return code %d",returnCode); |
---|
2903 | if ((returnCode&2) != 0) { |
---|
2904 | // could add cut |
---|
2905 | returnCode &= ~2; |
---|
2906 | //printf("could add cut with %d elements (if all 0-1)\n",nFix); |
---|
2907 | } else { |
---|
2908 | //printf("\n"); |
---|
2909 | } |
---|
2910 | } |
---|
2911 | fixPriority_ = -1; // switch off |
---|
2912 | |
---|
2913 | delete newSolver; |
---|
2914 | return returnCode; |
---|
2915 | } |
---|
2916 | // update model |
---|
2917 | void CbcHeuristicPartial::setModel(CbcModel * model) |
---|
2918 | { |
---|
2919 | model_ = model; |
---|
2920 | assert(model_->solver()); |
---|
2921 | // make sure model okay for heuristic |
---|
2922 | validate(); |
---|
2923 | } |
---|
2924 | // Validate model i.e. sets when_ to 0 if necessary (may be NULL) |
---|
2925 | void |
---|
2926 | CbcHeuristicPartial::validate() |
---|
2927 | { |
---|
2928 | if (model_ && (when() % 100) < 10) { |
---|
2929 | if (model_->numberIntegers() != |
---|
2930 | model_->numberObjects()) |
---|
2931 | setWhen(0); |
---|
2932 | } |
---|
2933 | } |
---|
2934 | bool |
---|
2935 | CbcHeuristicPartial::shouldHeurRun(int /*whereFrom*/) |
---|
2936 | { |
---|
2937 | return true; |
---|
2938 | } |
---|
2939 | |
---|
2940 | // Default Constructor |
---|
2941 | CbcSerendipity::CbcSerendipity() |
---|
2942 | : CbcHeuristic() |
---|
2943 | { |
---|
2944 | } |
---|
2945 | |
---|
2946 | // Constructor from model |
---|
2947 | CbcSerendipity::CbcSerendipity(CbcModel & model) |
---|
2948 | : CbcHeuristic(model) |
---|
2949 | { |
---|
2950 | } |
---|
2951 | |
---|
2952 | // Destructor |
---|
2953 | CbcSerendipity::~CbcSerendipity () |
---|
2954 | { |
---|
2955 | } |
---|
2956 | |
---|
2957 | // Clone |
---|
2958 | CbcHeuristic * |
---|
2959 | CbcSerendipity::clone() const |
---|
2960 | { |
---|
2961 | return new CbcSerendipity(*this); |
---|
2962 | } |
---|
2963 | // Create C++ lines to get to current state |
---|
2964 | void |
---|
2965 | CbcSerendipity::generateCpp( FILE * fp) |
---|
2966 | { |
---|
2967 | fprintf(fp, "0#include \"CbcHeuristic.hpp\"\n"); |
---|
2968 | fprintf(fp, "3 CbcSerendipity serendipity(*cbcModel);\n"); |
---|
2969 | CbcHeuristic::generateCpp(fp, "serendipity"); |
---|
2970 | fprintf(fp, "3 cbcModel->addHeuristic(&serendipity);\n"); |
---|
2971 | } |
---|
2972 | |
---|
2973 | // Copy constructor |
---|
2974 | CbcSerendipity::CbcSerendipity(const CbcSerendipity & rhs) |
---|
2975 | : |
---|
2976 | CbcHeuristic(rhs) |
---|
2977 | { |
---|
2978 | } |
---|
2979 | |
---|
2980 | // Assignment operator |
---|
2981 | CbcSerendipity & |
---|
2982 | CbcSerendipity::operator=( const CbcSerendipity & rhs) |
---|
2983 | { |
---|
2984 | if (this != &rhs) { |
---|
2985 | CbcHeuristic::operator=(rhs); |
---|
2986 | } |
---|
2987 | return *this; |
---|
2988 | } |
---|
2989 | |
---|
2990 | // Returns 1 if solution, 0 if not |
---|
2991 | int |
---|
2992 | CbcSerendipity::solution(double & solutionValue, |
---|
2993 | double * betterSolution) |
---|
2994 | { |
---|
2995 | if (!model_) |
---|
2996 | return 0; |
---|
2997 | #ifdef HEURISTIC_INFORM |
---|
2998 | printf("Entering heuristic %s - nRuns %d numCould %d when %d\n", |
---|
2999 | heuristicName(),numRuns_,numCouldRun_,when_); |
---|
3000 | #endif |
---|
3001 | if (!inputSolution_) { |
---|
3002 | // get information on solver type |
---|
3003 | OsiAuxInfo * auxInfo = model_->solver()->getAuxiliaryInfo(); |
---|
3004 | OsiBabSolver * auxiliaryInfo = dynamic_cast< OsiBabSolver *> (auxInfo); |
---|
3005 | if (auxiliaryInfo) { |
---|
3006 | return auxiliaryInfo->solution(solutionValue, betterSolution, model_->solver()->getNumCols()); |
---|
3007 | } else { |
---|
3008 | return 0; |
---|
3009 | } |
---|
3010 | } else { |
---|
3011 | int numberColumns = model_->getNumCols(); |
---|
3012 | double value = inputSolution_[numberColumns]; |
---|
3013 | int returnCode = 0; |
---|
3014 | if (value < solutionValue) { |
---|
3015 | solutionValue = value; |
---|
3016 | memcpy(betterSolution, inputSolution_, numberColumns*sizeof(double)); |
---|
3017 | returnCode = 1; |
---|
3018 | } |
---|
3019 | delete [] inputSolution_; |
---|
3020 | inputSolution_ = NULL; |
---|
3021 | model_ = NULL; // switch off |
---|
3022 | return returnCode; |
---|
3023 | } |
---|
3024 | } |
---|
3025 | // update model |
---|
3026 | void CbcSerendipity::setModel(CbcModel * model) |
---|
3027 | { |
---|
3028 | model_ = model; |
---|
3029 | } |
---|
3030 | // Resets stuff if model changes |
---|
3031 | void |
---|
3032 | CbcSerendipity::resetModel(CbcModel * model) |
---|
3033 | { |
---|
3034 | model_ = model; |
---|
3035 | } |
---|
3036 | |
---|
3037 | |
---|
3038 | // Default Constructor |
---|
3039 | CbcHeuristicJustOne::CbcHeuristicJustOne() |
---|
3040 | : CbcHeuristic(), |
---|
3041 | probabilities_(NULL), |
---|
3042 | heuristic_(NULL), |
---|
3043 | numberHeuristics_(0) |
---|
3044 | { |
---|
3045 | } |
---|
3046 | |
---|
3047 | // Constructor from model |
---|
3048 | CbcHeuristicJustOne::CbcHeuristicJustOne(CbcModel & model) |
---|
3049 | : CbcHeuristic(model), |
---|
3050 | probabilities_(NULL), |
---|
3051 | heuristic_(NULL), |
---|
3052 | numberHeuristics_(0) |
---|
3053 | { |
---|
3054 | } |
---|
3055 | |
---|
3056 | // Destructor |
---|
3057 | CbcHeuristicJustOne::~CbcHeuristicJustOne () |
---|
3058 | { |
---|
3059 | for (int i = 0; i < numberHeuristics_; i++) |
---|
3060 | delete heuristic_[i]; |
---|
3061 | delete [] heuristic_; |
---|
3062 | delete [] probabilities_; |
---|
3063 | } |
---|
3064 | |
---|
3065 | // Clone |
---|
3066 | CbcHeuristicJustOne * |
---|
3067 | CbcHeuristicJustOne::clone() const |
---|
3068 | { |
---|
3069 | return new CbcHeuristicJustOne(*this); |
---|
3070 | } |
---|
3071 | |
---|
3072 | // Create C++ lines to get to current state |
---|
3073 | void |
---|
3074 | CbcHeuristicJustOne::generateCpp( FILE * fp) |
---|
3075 | { |
---|
3076 | CbcHeuristicJustOne other; |
---|
3077 | fprintf(fp, "0#include \"CbcHeuristicJustOne.hpp\"\n"); |
---|
3078 | fprintf(fp, "3 CbcHeuristicJustOne heuristicJustOne(*cbcModel);\n"); |
---|
3079 | CbcHeuristic::generateCpp(fp, "heuristicJustOne"); |
---|
3080 | fprintf(fp, "3 cbcModel->addHeuristic(&heuristicJustOne);\n"); |
---|
3081 | } |
---|
3082 | |
---|
3083 | // Copy constructor |
---|
3084 | CbcHeuristicJustOne::CbcHeuristicJustOne(const CbcHeuristicJustOne & rhs) |
---|
3085 | : |
---|
3086 | CbcHeuristic(rhs), |
---|
3087 | probabilities_(NULL), |
---|
3088 | heuristic_(NULL), |
---|
3089 | numberHeuristics_(rhs.numberHeuristics_) |
---|
3090 | { |
---|
3091 | if (numberHeuristics_) { |
---|
3092 | probabilities_ = CoinCopyOfArray(rhs.probabilities_, numberHeuristics_); |
---|
3093 | heuristic_ = new CbcHeuristic * [numberHeuristics_]; |
---|
3094 | for (int i = 0; i < numberHeuristics_; i++) |
---|
3095 | heuristic_[i] = rhs.heuristic_[i]->clone(); |
---|
3096 | } |
---|
3097 | } |
---|
3098 | |
---|
3099 | // Assignment operator |
---|
3100 | CbcHeuristicJustOne & |
---|
3101 | CbcHeuristicJustOne::operator=( const CbcHeuristicJustOne & rhs) |
---|
3102 | { |
---|
3103 | if (this != &rhs) { |
---|
3104 | CbcHeuristic::operator=(rhs); |
---|
3105 | for (int i = 0; i < numberHeuristics_; i++) |
---|
3106 | delete heuristic_[i]; |
---|
3107 | delete [] heuristic_; |
---|
3108 | delete [] probabilities_; |
---|
3109 | probabilities_ = NULL; |
---|
3110 | heuristic_ = NULL; |
---|
3111 | numberHeuristics_ = rhs.numberHeuristics_; |
---|
3112 | if (numberHeuristics_) { |
---|
3113 | probabilities_ = CoinCopyOfArray(rhs.probabilities_, numberHeuristics_); |
---|
3114 | heuristic_ = new CbcHeuristic * [numberHeuristics_]; |
---|
3115 | for (int i = 0; i < numberHeuristics_; i++) |
---|
3116 | heuristic_[i] = rhs.heuristic_[i]->clone(); |
---|
3117 | } |
---|
3118 | } |
---|
3119 | return *this; |
---|
3120 | } |
---|
3121 | // Sets value of solution |
---|
3122 | // Returns 1 if solution, 0 if not |
---|
3123 | int |
---|
3124 | CbcHeuristicJustOne::solution(double & solutionValue, |
---|
3125 | double * betterSolution) |
---|
3126 | { |
---|
3127 | #ifdef DIVE_DEBUG |
---|
3128 | std::cout << "solutionValue = " << solutionValue << std::endl; |
---|
3129 | #endif |
---|
3130 | ++numCouldRun_; |
---|
3131 | |
---|
3132 | // test if the heuristic can run |
---|
3133 | if (!shouldHeurRun_randomChoice() || !numberHeuristics_) |
---|
3134 | return 0; |
---|
3135 | double randomNumber = randomNumberGenerator_.randomDouble(); |
---|
3136 | int i; |
---|
3137 | for (i = 0; i < numberHeuristics_; i++) { |
---|
3138 | if (randomNumber < probabilities_[i]) |
---|
3139 | break; |
---|
3140 | } |
---|
3141 | assert (i < numberHeuristics_); |
---|
3142 | int returnCode; |
---|
3143 | //model_->unsetDivingHasRun(); |
---|
3144 | #ifdef COIN_DEVELOP |
---|
3145 | printf("JustOne running %s\n", |
---|
3146 | heuristic_[i]->heuristicName()); |
---|
3147 | #endif |
---|
3148 | returnCode = heuristic_[i]->solution(solutionValue, betterSolution); |
---|
3149 | #ifdef COIN_DEVELOP |
---|
3150 | if (returnCode) |
---|
3151 | printf("JustOne running %s found solution\n", |
---|
3152 | heuristic_[i]->heuristicName()); |
---|
3153 | #endif |
---|
3154 | return returnCode; |
---|
3155 | } |
---|
3156 | // Resets stuff if model changes |
---|
3157 | void |
---|
3158 | CbcHeuristicJustOne::resetModel(CbcModel * model) |
---|
3159 | { |
---|
3160 | CbcHeuristic::resetModel(model); |
---|
3161 | for (int i = 0; i < numberHeuristics_; i++) |
---|
3162 | heuristic_[i]->resetModel(model); |
---|
3163 | } |
---|
3164 | // update model (This is needed if cliques update matrix etc) |
---|
3165 | void |
---|
3166 | CbcHeuristicJustOne::setModel(CbcModel * model) |
---|
3167 | { |
---|
3168 | CbcHeuristic::setModel(model); |
---|
3169 | for (int i = 0; i < numberHeuristics_; i++) |
---|
3170 | heuristic_[i]->setModel(model); |
---|
3171 | } |
---|
3172 | // Validate model i.e. sets when_ to 0 if necessary (may be NULL) |
---|
3173 | void |
---|
3174 | CbcHeuristicJustOne::validate() |
---|
3175 | { |
---|
3176 | CbcHeuristic::validate(); |
---|
3177 | for (int i = 0; i < numberHeuristics_; i++) |
---|
3178 | heuristic_[i]->validate(); |
---|
3179 | } |
---|
3180 | // Adds an heuristic with probability |
---|
3181 | void |
---|
3182 | CbcHeuristicJustOne::addHeuristic(const CbcHeuristic * heuristic, double probability) |
---|
3183 | { |
---|
3184 | CbcHeuristic * thisOne = heuristic->clone(); |
---|
3185 | thisOne->setWhen(-999); |
---|
3186 | CbcHeuristic ** tempH = CoinCopyOfArrayPartial(heuristic_, numberHeuristics_ + 1, |
---|
3187 | numberHeuristics_); |
---|
3188 | delete [] heuristic_; |
---|
3189 | heuristic_ = tempH; |
---|
3190 | heuristic_[numberHeuristics_] = thisOne; |
---|
3191 | double * tempP = CoinCopyOfArrayPartial(probabilities_, numberHeuristics_ + 1, |
---|
3192 | numberHeuristics_); |
---|
3193 | delete [] probabilities_; |
---|
3194 | probabilities_ = tempP; |
---|
3195 | probabilities_[numberHeuristics_] = probability; |
---|
3196 | numberHeuristics_++; |
---|
3197 | } |
---|
3198 | // Normalize probabilities |
---|
3199 | void |
---|
3200 | CbcHeuristicJustOne::normalizeProbabilities() |
---|
3201 | { |
---|
3202 | double sum = 0.0; |
---|
3203 | for (int i = 0; i < numberHeuristics_; i++) |
---|
3204 | sum += probabilities_[i]; |
---|
3205 | double multiplier = 1.0 / sum; |
---|
3206 | sum = 0.0; |
---|
3207 | for (int i = 0; i < numberHeuristics_; i++) { |
---|
3208 | sum += probabilities_[i]; |
---|
3209 | probabilities_[i] = sum * multiplier; |
---|
3210 | } |
---|
3211 | assert (fabs(probabilities_[numberHeuristics_-1] - 1.0) < 1.0e-5); |
---|
3212 | probabilities_[numberHeuristics_-1] = 1.000001; |
---|
3213 | } |
---|
3214 | |
---|