1 | /* $Id: CbcHeuristicLocal.cpp 2105 2015-01-05 13:11:11Z 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 | #include <cassert> |
---|
11 | #include <cstdlib> |
---|
12 | #include <cmath> |
---|
13 | #include <cfloat> |
---|
14 | |
---|
15 | #include "OsiSolverInterface.hpp" |
---|
16 | #include "CbcModel.hpp" |
---|
17 | #include "CbcMessage.hpp" |
---|
18 | #include "CbcHeuristicLocal.hpp" |
---|
19 | #include "CbcHeuristicFPump.hpp" |
---|
20 | #include "CbcBranchActual.hpp" |
---|
21 | #include "CbcStrategy.hpp" |
---|
22 | #include "CglPreProcess.hpp" |
---|
23 | |
---|
24 | // Default Constructor |
---|
25 | CbcHeuristicLocal::CbcHeuristicLocal() |
---|
26 | : CbcHeuristic() |
---|
27 | { |
---|
28 | numberSolutions_ = 0; |
---|
29 | swap_ = 0; |
---|
30 | used_ = NULL; |
---|
31 | lastRunDeep_ = -1000000; |
---|
32 | switches_ |= 16; // needs a new solution |
---|
33 | } |
---|
34 | |
---|
35 | // Constructor with model - assumed before cuts |
---|
36 | |
---|
37 | CbcHeuristicLocal::CbcHeuristicLocal(CbcModel & model) |
---|
38 | : CbcHeuristic(model) |
---|
39 | { |
---|
40 | numberSolutions_ = 0; |
---|
41 | swap_ = 0; |
---|
42 | lastRunDeep_ = -1000000; |
---|
43 | switches_ |= 16; // needs a new solution |
---|
44 | // Get a copy of original matrix |
---|
45 | assert(model.solver()); |
---|
46 | if (model.solver()->getNumRows()) { |
---|
47 | matrix_ = *model.solver()->getMatrixByCol(); |
---|
48 | } |
---|
49 | int numberColumns = model.solver()->getNumCols(); |
---|
50 | used_ = new int[numberColumns]; |
---|
51 | memset(used_, 0, numberColumns*sizeof(int)); |
---|
52 | } |
---|
53 | |
---|
54 | // Destructor |
---|
55 | CbcHeuristicLocal::~CbcHeuristicLocal () |
---|
56 | { |
---|
57 | delete [] used_; |
---|
58 | } |
---|
59 | |
---|
60 | // Clone |
---|
61 | CbcHeuristic * |
---|
62 | CbcHeuristicLocal::clone() const |
---|
63 | { |
---|
64 | return new CbcHeuristicLocal(*this); |
---|
65 | } |
---|
66 | // Create C++ lines to get to current state |
---|
67 | void |
---|
68 | CbcHeuristicLocal::generateCpp( FILE * fp) |
---|
69 | { |
---|
70 | CbcHeuristicLocal other; |
---|
71 | fprintf(fp, "0#include \"CbcHeuristicLocal.hpp\"\n"); |
---|
72 | fprintf(fp, "3 CbcHeuristicLocal heuristicLocal(*cbcModel);\n"); |
---|
73 | CbcHeuristic::generateCpp(fp, "heuristicLocal"); |
---|
74 | if (swap_ != other.swap_) |
---|
75 | fprintf(fp, "3 heuristicLocal.setSearchType(%d);\n", swap_); |
---|
76 | else |
---|
77 | fprintf(fp, "4 heuristicLocal.setSearchType(%d);\n", swap_); |
---|
78 | fprintf(fp, "3 cbcModel->addHeuristic(&heuristicLocal);\n"); |
---|
79 | } |
---|
80 | |
---|
81 | // Copy constructor |
---|
82 | CbcHeuristicLocal::CbcHeuristicLocal(const CbcHeuristicLocal & rhs) |
---|
83 | : |
---|
84 | CbcHeuristic(rhs), |
---|
85 | matrix_(rhs.matrix_), |
---|
86 | numberSolutions_(rhs.numberSolutions_), |
---|
87 | swap_(rhs.swap_) |
---|
88 | { |
---|
89 | if (model_ && rhs.used_) { |
---|
90 | int numberColumns = model_->solver()->getNumCols(); |
---|
91 | used_ = CoinCopyOfArray(rhs.used_, numberColumns); |
---|
92 | } else { |
---|
93 | used_ = NULL; |
---|
94 | } |
---|
95 | } |
---|
96 | |
---|
97 | // Assignment operator |
---|
98 | CbcHeuristicLocal & |
---|
99 | CbcHeuristicLocal::operator=( const CbcHeuristicLocal & rhs) |
---|
100 | { |
---|
101 | if (this != &rhs) { |
---|
102 | CbcHeuristic::operator=(rhs); |
---|
103 | matrix_ = rhs.matrix_; |
---|
104 | numberSolutions_ = rhs.numberSolutions_; |
---|
105 | swap_ = rhs.swap_; |
---|
106 | delete [] used_; |
---|
107 | if (model_ && rhs.used_) { |
---|
108 | int numberColumns = model_->solver()->getNumCols(); |
---|
109 | used_ = CoinCopyOfArray(rhs.used_, numberColumns); |
---|
110 | } else { |
---|
111 | used_ = NULL; |
---|
112 | } |
---|
113 | } |
---|
114 | return *this; |
---|
115 | } |
---|
116 | |
---|
117 | // Resets stuff if model changes |
---|
118 | void |
---|
119 | CbcHeuristicLocal::resetModel(CbcModel * /*model*/) |
---|
120 | { |
---|
121 | //CbcHeuristic::resetModel(model); |
---|
122 | delete [] used_; |
---|
123 | if (model_ && used_) { |
---|
124 | int numberColumns = model_->solver()->getNumCols(); |
---|
125 | used_ = new int[numberColumns]; |
---|
126 | memset(used_, 0, numberColumns*sizeof(int)); |
---|
127 | } else { |
---|
128 | used_ = NULL; |
---|
129 | } |
---|
130 | } |
---|
131 | /* |
---|
132 | Run a mini-BaB search after fixing all variables not marked as used by |
---|
133 | solution(). (See comments there for semantics.) |
---|
134 | |
---|
135 | Return values are: |
---|
136 | 1: smallBranchAndBound found a solution |
---|
137 | 0: everything else |
---|
138 | |
---|
139 | The degree of overload as return codes from smallBranchAndBound are folded |
---|
140 | into 0 is such that it's impossible to distinguish return codes that really |
---|
141 | require attention from a simple `nothing of interest'. |
---|
142 | */ |
---|
143 | // This version fixes stuff and does IP |
---|
144 | int |
---|
145 | CbcHeuristicLocal::solutionFix(double & objectiveValue, |
---|
146 | double * newSolution, |
---|
147 | const int * /*keep*/) |
---|
148 | { |
---|
149 | /* |
---|
150 | If when is set to off (0), or set to root (1) and we're not at the root, |
---|
151 | return. If this heuristic discovered the current solution, don't continue. |
---|
152 | */ |
---|
153 | |
---|
154 | numCouldRun_++; |
---|
155 | // See if to do |
---|
156 | if (!when() || (when() == 1 && model_->phase() != 1)) |
---|
157 | return 0; // switched off |
---|
158 | // Don't do if it was this heuristic which found solution! |
---|
159 | if (this == model_->lastHeuristic()) |
---|
160 | return 0; |
---|
161 | /* |
---|
162 | Load up a new solver with the solution. |
---|
163 | |
---|
164 | Why continuousSolver(), as opposed to solver()? |
---|
165 | */ |
---|
166 | OsiSolverInterface * newSolver = model_->continuousSolver()->clone(); |
---|
167 | const double * colLower = newSolver->getColLower(); |
---|
168 | //const double * colUpper = newSolver->getColUpper(); |
---|
169 | |
---|
170 | int numberIntegers = model_->numberIntegers(); |
---|
171 | const int * integerVariable = model_->integerVariable(); |
---|
172 | /* |
---|
173 | The net effect here is that anything that hasn't moved from its lower bound |
---|
174 | will be fixed at lower bound. |
---|
175 | |
---|
176 | See comments in solution() w.r.t. asymmetric treatment of upper and lower |
---|
177 | bounds. |
---|
178 | */ |
---|
179 | |
---|
180 | int i; |
---|
181 | int nFix = 0; |
---|
182 | for (i = 0; i < numberIntegers; i++) { |
---|
183 | int iColumn = integerVariable[i]; |
---|
184 | const OsiObject * object = model_->object(i); |
---|
185 | // get original bounds |
---|
186 | double originalLower; |
---|
187 | double originalUpper; |
---|
188 | getIntegerInformation( object, originalLower, originalUpper); |
---|
189 | newSolver->setColLower(iColumn, CoinMax(colLower[iColumn], originalLower)); |
---|
190 | if (!used_[iColumn]) { |
---|
191 | newSolver->setColUpper(iColumn, colLower[iColumn]); |
---|
192 | nFix++; |
---|
193 | } |
---|
194 | } |
---|
195 | /* |
---|
196 | Try a `small' branch-and-bound search. The notion here is that we've fixed a |
---|
197 | lot of variables and reduced the amount of `free' problem to a point where a |
---|
198 | small BaB search will suffice to fully explore the remaining problem. This |
---|
199 | routine will execute integer presolve, then call branchAndBound to do the |
---|
200 | actual search. |
---|
201 | */ |
---|
202 | int returnCode = 0; |
---|
203 | #ifdef CLP_INVESTIGATE2 |
---|
204 | printf("Fixing %d out of %d (%d continuous)\n", |
---|
205 | nFix, numberIntegers, newSolver->getNumCols() - numberIntegers); |
---|
206 | #endif |
---|
207 | if (nFix*10 <= numberIntegers) { |
---|
208 | // see if we can fix more |
---|
209 | int * which = new int [2*(numberIntegers-nFix)]; |
---|
210 | int * sort = which + (numberIntegers - nFix); |
---|
211 | int n = 0; |
---|
212 | for (i = 0; i < numberIntegers; i++) { |
---|
213 | int iColumn = integerVariable[i]; |
---|
214 | if (used_[iColumn]) { |
---|
215 | which[n] = iColumn; |
---|
216 | sort[n++] = used_[iColumn]; |
---|
217 | } |
---|
218 | } |
---|
219 | CoinSort_2(sort, sort + n, which); |
---|
220 | // only half fixed in total |
---|
221 | n = CoinMin(n, numberIntegers / 2 - nFix); |
---|
222 | int allow = CoinMax(numberSolutions_ - 2, sort[0]); |
---|
223 | int nFix2 = 0; |
---|
224 | for (i = 0; i < n; i++) { |
---|
225 | int iColumn = integerVariable[i]; |
---|
226 | if (used_[iColumn] <= allow) { |
---|
227 | newSolver->setColUpper(iColumn, colLower[iColumn]); |
---|
228 | nFix2++; |
---|
229 | } else { |
---|
230 | break; |
---|
231 | } |
---|
232 | } |
---|
233 | delete [] which; |
---|
234 | nFix += nFix2; |
---|
235 | #ifdef CLP_INVESTIGATE2 |
---|
236 | printf("Number fixed increased from %d to %d\n", |
---|
237 | nFix - nFix2, nFix); |
---|
238 | #endif |
---|
239 | } |
---|
240 | if (nFix*10 > numberIntegers) { |
---|
241 | returnCode = smallBranchAndBound(newSolver, numberNodes_, newSolution, objectiveValue, |
---|
242 | objectiveValue, "CbcHeuristicLocal"); |
---|
243 | /* |
---|
244 | -2 is return due to user event, and -1 is overloaded with what look to be |
---|
245 | two contradictory meanings. |
---|
246 | */ |
---|
247 | if (returnCode < 0) { |
---|
248 | returnCode = 0; // returned on size |
---|
249 | int numberColumns = newSolver->getNumCols(); |
---|
250 | int numberContinuous = numberColumns - numberIntegers; |
---|
251 | if (numberContinuous > 2*numberIntegers && |
---|
252 | nFix*10 < numberColumns) { |
---|
253 | #define LOCAL_FIX_CONTINUOUS |
---|
254 | #ifdef LOCAL_FIX_CONTINUOUS |
---|
255 | //const double * colUpper = newSolver->getColUpper(); |
---|
256 | const double * colLower = newSolver->getColLower(); |
---|
257 | int nAtLb = 0; |
---|
258 | //double sumDj=0.0; |
---|
259 | const double * dj = newSolver->getReducedCost(); |
---|
260 | double direction = newSolver->getObjSense(); |
---|
261 | for (int iColumn = 0; iColumn < numberColumns; iColumn++) { |
---|
262 | if (!newSolver->isInteger(iColumn)) { |
---|
263 | if (!used_[iColumn]) { |
---|
264 | //double djValue = dj[iColumn]*direction; |
---|
265 | nAtLb++; |
---|
266 | //sumDj += djValue; |
---|
267 | } |
---|
268 | } |
---|
269 | } |
---|
270 | if (nAtLb) { |
---|
271 | // fix some continuous |
---|
272 | double * sort = new double[nAtLb]; |
---|
273 | int * which = new int [nAtLb]; |
---|
274 | //double threshold = CoinMax((0.01*sumDj)/static_cast<double>(nAtLb),1.0e-6); |
---|
275 | int nFix2 = 0; |
---|
276 | for (int iColumn = 0; iColumn < numberColumns; iColumn++) { |
---|
277 | if (!newSolver->isInteger(iColumn)) { |
---|
278 | if (!used_[iColumn]) { |
---|
279 | double djValue = dj[iColumn] * direction; |
---|
280 | if (djValue > 1.0e-6) { |
---|
281 | sort[nFix2] = -djValue; |
---|
282 | which[nFix2++] = iColumn; |
---|
283 | } |
---|
284 | } |
---|
285 | } |
---|
286 | } |
---|
287 | CoinSort_2(sort, sort + nFix2, which); |
---|
288 | int divisor = 2; |
---|
289 | nFix2 = CoinMin(nFix2, (numberColumns - nFix) / divisor); |
---|
290 | for (int i = 0; i < nFix2; i++) { |
---|
291 | int iColumn = which[i]; |
---|
292 | newSolver->setColUpper(iColumn, colLower[iColumn]); |
---|
293 | } |
---|
294 | delete [] sort; |
---|
295 | delete [] which; |
---|
296 | #ifdef CLP_INVESTIGATE2 |
---|
297 | printf("%d integers have zero value, and %d continuous fixed at lb\n", |
---|
298 | nFix, nFix2); |
---|
299 | #endif |
---|
300 | returnCode = smallBranchAndBound(newSolver, |
---|
301 | numberNodes_, newSolution, |
---|
302 | objectiveValue, |
---|
303 | objectiveValue, "CbcHeuristicLocal"); |
---|
304 | if (returnCode < 0) |
---|
305 | returnCode = 0; // returned on size |
---|
306 | } |
---|
307 | #endif |
---|
308 | } |
---|
309 | } |
---|
310 | } |
---|
311 | /* |
---|
312 | If the result is complete exploration with a solution (3) or proven |
---|
313 | infeasibility (2), we could generate a cut (the AI folks would call it a |
---|
314 | nogood) to prevent us from going down this route in the future. |
---|
315 | */ |
---|
316 | if ((returnCode&2) != 0) { |
---|
317 | // could add cut |
---|
318 | returnCode &= ~2; |
---|
319 | } |
---|
320 | |
---|
321 | delete newSolver; |
---|
322 | return returnCode; |
---|
323 | } |
---|
324 | /* |
---|
325 | First tries setting a variable to better value. If feasible then |
---|
326 | tries setting others. If not feasible then tries swaps |
---|
327 | Returns 1 if solution, 0 if not |
---|
328 | The main body of this routine implements an O((q^2)/2) brute force search |
---|
329 | around the current solution, for q = number of integer variables. Call this |
---|
330 | the inc/dec heuristic. For each integer variable x<i>, first decrement the |
---|
331 | value. Then, for integer variables x<i+1>, ..., x<q-1>, try increment and |
---|
332 | decrement. If one of these permutations produces a better solution, |
---|
333 | remember it. Then repeat, with x<i> incremented. If we find a better |
---|
334 | solution, update our notion of current solution and continue. |
---|
335 | |
---|
336 | The net effect is a greedy walk: As each improving pair is found, the |
---|
337 | current solution is updated and the search continues from this updated |
---|
338 | solution. |
---|
339 | |
---|
340 | Way down at the end, we call solutionFix, which will create a drastically |
---|
341 | restricted problem based on variables marked as used, then do mini-BaC on |
---|
342 | the restricted problem. This can occur even if we don't try the inc/dec |
---|
343 | heuristic. This would be more obvious if the inc/dec heuristic were broken |
---|
344 | out as a separate routine and solutionFix had a name that reflected where |
---|
345 | it was headed. |
---|
346 | |
---|
347 | The return code of 0 is grossly overloaded, because it maps to a return |
---|
348 | code of 0 from solutionFix, which is itself grossly overloaded. See |
---|
349 | comments in solutionFix and in CbcHeuristic::smallBranchAndBound. |
---|
350 | */ |
---|
351 | int |
---|
352 | CbcHeuristicLocal::solution(double & solutionValue, |
---|
353 | double * betterSolution) |
---|
354 | { |
---|
355 | /* |
---|
356 | Execute only if a new solution has been discovered since the last time we |
---|
357 | were called. |
---|
358 | */ |
---|
359 | |
---|
360 | numCouldRun_++; |
---|
361 | // See if frequency kills off idea |
---|
362 | int swap = swap_%100; |
---|
363 | int skip = swap_/100; |
---|
364 | int nodeCount = model_->getNodeCount(); |
---|
365 | if (nodeCount<lastRunDeep_+skip && nodeCount != lastRunDeep_+1) |
---|
366 | return 0; |
---|
367 | if (numberSolutions_ == model_->getSolutionCount() && |
---|
368 | (numberSolutions_ == howOftenShallow_ || |
---|
369 | nodeCount < lastRunDeep_+2*skip)) |
---|
370 | return 0; |
---|
371 | howOftenShallow_ = numberSolutions_; |
---|
372 | numberSolutions_ = model_->getSolutionCount(); |
---|
373 | if (nodeCount<lastRunDeep_+skip ) |
---|
374 | return 0; |
---|
375 | #ifdef HEURISTIC_INFORM |
---|
376 | printf("Entering heuristic %s - nRuns %d numCould %d when %d\n", |
---|
377 | heuristicName(),numRuns_,numCouldRun_,when_); |
---|
378 | #endif |
---|
379 | lastRunDeep_ = nodeCount; |
---|
380 | howOftenShallow_ = numberSolutions_; |
---|
381 | |
---|
382 | if ((swap%10) == 2) { |
---|
383 | // try merge |
---|
384 | return solutionFix( solutionValue, betterSolution, NULL); |
---|
385 | } |
---|
386 | /* |
---|
387 | Exclude long (column), thin (row) systems. |
---|
388 | |
---|
389 | Given the n^2 nature of the search, more than 100,000 columns could get |
---|
390 | expensive. But I don't yet see the rationale for the second part of the |
---|
391 | condition (cols > 10*rows). And cost is proportional to number of integer |
---|
392 | variables --- shouldn't we use that? |
---|
393 | |
---|
394 | Why wait until we have more than one solution? |
---|
395 | */ |
---|
396 | if ((model_->getNumCols() > 100000 && model_->getNumCols() > |
---|
397 | 10*model_->getNumRows()) || numberSolutions_ <= 1) |
---|
398 | return 0; // probably not worth it |
---|
399 | // worth trying |
---|
400 | |
---|
401 | OsiSolverInterface * solver = model_->solver(); |
---|
402 | const double * rowLower = solver->getRowLower(); |
---|
403 | const double * rowUpper = solver->getRowUpper(); |
---|
404 | const double * solution = model_->bestSolution(); |
---|
405 | /* |
---|
406 | Shouldn't this test be redundant if we've already checked that |
---|
407 | numberSolutions_ > 1? Stronger: shouldn't this be an assertion? |
---|
408 | */ |
---|
409 | if (!solution) |
---|
410 | return 0; // No solution found yet |
---|
411 | const double * objective = solver->getObjCoefficients(); |
---|
412 | double primalTolerance; |
---|
413 | solver->getDblParam(OsiPrimalTolerance, primalTolerance); |
---|
414 | |
---|
415 | int numberRows = matrix_.getNumRows(); |
---|
416 | |
---|
417 | int numberIntegers = model_->numberIntegers(); |
---|
418 | const int * integerVariable = model_->integerVariable(); |
---|
419 | |
---|
420 | int i; |
---|
421 | double direction = solver->getObjSense(); |
---|
422 | double newSolutionValue = model_->getObjValue() * direction; |
---|
423 | int returnCode = 0; |
---|
424 | numRuns_++; |
---|
425 | // Column copy |
---|
426 | const double * element = matrix_.getElements(); |
---|
427 | const int * row = matrix_.getIndices(); |
---|
428 | const CoinBigIndex * columnStart = matrix_.getVectorStarts(); |
---|
429 | const int * columnLength = matrix_.getVectorLengths(); |
---|
430 | |
---|
431 | // Get solution array for heuristic solution |
---|
432 | int numberColumns = solver->getNumCols(); |
---|
433 | double * newSolution = new double [numberColumns]; |
---|
434 | memcpy(newSolution, solution, numberColumns*sizeof(double)); |
---|
435 | #ifdef LOCAL_FIX_CONTINUOUS |
---|
436 | // mark continuous used |
---|
437 | const double * columnLower = solver->getColLower(); |
---|
438 | for (int iColumn = 0; iColumn < numberColumns; iColumn++) { |
---|
439 | if (!solver->isInteger(iColumn)) { |
---|
440 | if (solution[iColumn] > columnLower[iColumn] + 1.0e-8) |
---|
441 | used_[iColumn] = numberSolutions_; |
---|
442 | } |
---|
443 | } |
---|
444 | #endif |
---|
445 | |
---|
446 | // way is 1 if down possible, 2 if up possible, 3 if both possible |
---|
447 | char * way = new char[numberIntegers]; |
---|
448 | // corrected costs |
---|
449 | double * cost = new double[numberIntegers]; |
---|
450 | // for array to mark infeasible rows after iColumn branch |
---|
451 | char * mark = new char[numberRows]; |
---|
452 | memset(mark, 0, numberRows); |
---|
453 | // space to save values so we don't introduce rounding errors |
---|
454 | double * save = new double[numberRows]; |
---|
455 | /* |
---|
456 | Force variables within their original bounds, then to the nearest integer. |
---|
457 | Overall, we seem to be prepared to cope with noninteger bounds. Is this |
---|
458 | necessary? Seems like we'd be better off to force the bounds to integrality |
---|
459 | as part of preprocessing. More generally, why do we need to do this? This |
---|
460 | solution should have been cleaned and checked when it was accepted as a |
---|
461 | solution! |
---|
462 | |
---|
463 | Once the value is set, decide whether we can move up or down. |
---|
464 | |
---|
465 | The only place that used_ is used is in solutionFix; if a variable is not |
---|
466 | flagged as used, it will be fixed (at lower bound). Why the asymmetric |
---|
467 | treatment? This makes some sense for binary variables (for which there are |
---|
468 | only two options). But for general integer variables, why not make a similar |
---|
469 | test against the original upper bound? |
---|
470 | */ |
---|
471 | |
---|
472 | // clean solution |
---|
473 | for (i = 0; i < numberIntegers; i++) { |
---|
474 | int iColumn = integerVariable[i]; |
---|
475 | const OsiObject * object = model_->object(i); |
---|
476 | // get original bounds |
---|
477 | double originalLower; |
---|
478 | double originalUpper; |
---|
479 | getIntegerInformation( object, originalLower, originalUpper); |
---|
480 | double value = newSolution[iColumn]; |
---|
481 | if (value < originalLower) { |
---|
482 | value = originalLower; |
---|
483 | newSolution[iColumn] = value; |
---|
484 | } else if (value > originalUpper) { |
---|
485 | value = originalUpper; |
---|
486 | newSolution[iColumn] = value; |
---|
487 | } |
---|
488 | double nearest = floor(value + 0.5); |
---|
489 | //assert(fabs(value-nearest)<10.0*primalTolerance); |
---|
490 | value = nearest; |
---|
491 | newSolution[iColumn] = nearest; |
---|
492 | // if away from lower bound mark that fact |
---|
493 | if (nearest > originalLower) { |
---|
494 | used_[iColumn] = numberSolutions_; |
---|
495 | } |
---|
496 | cost[i] = direction * objective[iColumn]; |
---|
497 | /* |
---|
498 | Given previous computation we're checking that value is at least 1 away |
---|
499 | from the original bounds. |
---|
500 | */ |
---|
501 | int iway = 0; |
---|
502 | |
---|
503 | if (value > originalLower + 0.5) |
---|
504 | iway = 1; |
---|
505 | if (value < originalUpper - 0.5) |
---|
506 | iway |= 2; |
---|
507 | way[i] = static_cast<char>(iway); |
---|
508 | } |
---|
509 | /* |
---|
510 | Calculate lhs of each constraint for groomed solution. |
---|
511 | */ |
---|
512 | // get row activities |
---|
513 | double * rowActivity = new double[numberRows]; |
---|
514 | memset(rowActivity, 0, numberRows*sizeof(double)); |
---|
515 | |
---|
516 | for (i = 0; i < numberColumns; i++) { |
---|
517 | int j; |
---|
518 | double value = newSolution[i]; |
---|
519 | if (value) { |
---|
520 | for (j = columnStart[i]; |
---|
521 | j < columnStart[i] + columnLength[i]; j++) { |
---|
522 | int iRow = row[j]; |
---|
523 | rowActivity[iRow] += value * element[j]; |
---|
524 | } |
---|
525 | } |
---|
526 | } |
---|
527 | /* |
---|
528 | Check that constraints are satisfied. For small infeasibility, force the |
---|
529 | activity within bound. Again, why is this necessary if the current solution |
---|
530 | was accepted as a valid solution? |
---|
531 | |
---|
532 | Why are we scanning past the first unacceptable constraint? |
---|
533 | */ |
---|
534 | // check was feasible - if not adjust (cleaning may move) |
---|
535 | // if very infeasible then give up |
---|
536 | bool tryHeuristic = true; |
---|
537 | for (i = 0; i < numberRows; i++) { |
---|
538 | if (rowActivity[i] < rowLower[i]) { |
---|
539 | if (rowActivity[i] < rowLower[i] - 10.0*primalTolerance) |
---|
540 | tryHeuristic = false; |
---|
541 | rowActivity[i] = rowLower[i]; |
---|
542 | } else if (rowActivity[i] > rowUpper[i]) { |
---|
543 | if (rowActivity[i] < rowUpper[i] + 10.0*primalTolerance) |
---|
544 | tryHeuristic = false; |
---|
545 | rowActivity[i] = rowUpper[i]; |
---|
546 | } |
---|
547 | } |
---|
548 | /* |
---|
549 | This bit of code is not quite totally redundant: it'll bail at 10,000 |
---|
550 | instead of 100,000. Potentially we can do a lot of work to get here, only |
---|
551 | to abandon it. |
---|
552 | */ |
---|
553 | // Switch off if may take too long |
---|
554 | if (model_->getNumCols() > 10000 && model_->getNumCols() > |
---|
555 | 10*model_->getNumRows()&&swap<10) |
---|
556 | tryHeuristic = false; |
---|
557 | /* |
---|
558 | Try the inc/dec heuristic? |
---|
559 | */ |
---|
560 | if (tryHeuristic) { |
---|
561 | |
---|
562 | // total change in objective |
---|
563 | double totalChange = 0.0; |
---|
564 | // local best change in objective |
---|
565 | double bestChange = 0.0; |
---|
566 | // maybe just do 1000 |
---|
567 | int maxIntegers = numberIntegers; |
---|
568 | // stop if too many goes |
---|
569 | int maxTries=COIN_INT_MAX; |
---|
570 | // integerVariable may be randomized copy! |
---|
571 | int * integerVariable = |
---|
572 | CoinCopyOfArray(model_->integerVariable(),numberIntegers); |
---|
573 | if (swap>9 && numberIntegers>500) { |
---|
574 | int type=swap/10; |
---|
575 | if (type==1) { |
---|
576 | // reduce |
---|
577 | maxIntegers = CoinMin(1000,numberIntegers); |
---|
578 | } else if (type==2) { |
---|
579 | // reduce even more |
---|
580 | maxTries=100000; |
---|
581 | maxIntegers = CoinMin(500,numberIntegers); |
---|
582 | } else if (type>2) { |
---|
583 | assert (type<10); |
---|
584 | int totals[7]={1000,500,100,50,50,50,50}; |
---|
585 | maxIntegers=CoinMin(totals[type-3],numberIntegers); |
---|
586 | double * weight = new double[numberIntegers]; |
---|
587 | for (int i=0;i<numberIntegers;i++) { |
---|
588 | weight[i]=model_->randomNumberGenerator()->randomDouble(); |
---|
589 | } |
---|
590 | CoinSort_2(weight,weight+numberIntegers,integerVariable); |
---|
591 | } |
---|
592 | } |
---|
593 | /* |
---|
594 | Outer loop to walk integer variables. Call the current variable x<i>. At the |
---|
595 | end of this loop, bestChange will contain the best (negative) change in the |
---|
596 | objective for any single pair. |
---|
597 | |
---|
598 | The trouble is, we're limited to monotonically increasing improvement. |
---|
599 | Suppose we discover an improvement of 10 for some pair. If, later in the |
---|
600 | search, we discover an improvement of 9 for some other pair, we will not use |
---|
601 | it. That seems wasteful. |
---|
602 | */ |
---|
603 | |
---|
604 | for (i = 0; i < numberIntegers; i++) { |
---|
605 | int iColumn = integerVariable[i]; |
---|
606 | bestChange = 0.0; |
---|
607 | int endInner = CoinMin(numberIntegers,i+maxIntegers); |
---|
608 | |
---|
609 | double objectiveCoefficient = cost[i]; |
---|
610 | int k; |
---|
611 | int j; |
---|
612 | int goodK = -1; |
---|
613 | int wayK = -1, wayI = -1; |
---|
614 | /* |
---|
615 | Try decrementing x<i>. |
---|
616 | */ |
---|
617 | if ((way[i]&1) != 0) { |
---|
618 | int numberInfeasible = 0; |
---|
619 | /* |
---|
620 | Adjust row activities where x<i> has a nonzero coefficient. Save the old |
---|
621 | values for restoration. Mark any rows that become infeasible as a result |
---|
622 | of the decrement. |
---|
623 | */ |
---|
624 | // save row activities and adjust |
---|
625 | for (j = columnStart[iColumn]; |
---|
626 | j < columnStart[iColumn] + columnLength[iColumn]; j++) { |
---|
627 | int iRow = row[j]; |
---|
628 | save[iRow] = rowActivity[iRow]; |
---|
629 | rowActivity[iRow] -= element[j]; |
---|
630 | if (rowActivity[iRow] < rowLower[iRow] - primalTolerance || |
---|
631 | rowActivity[iRow] > rowUpper[iRow] + primalTolerance) { |
---|
632 | // mark row |
---|
633 | mark[iRow] = 1; |
---|
634 | numberInfeasible++; |
---|
635 | } |
---|
636 | } |
---|
637 | /* |
---|
638 | Run through the remaining integer variables. Try increment and decrement on |
---|
639 | each one. If the potential objective change is better than anything we've |
---|
640 | seen so far, do a full evaluation of x<k> in that direction. If we can |
---|
641 | repair all infeasibilities introduced by pushing x<i> down, we have a |
---|
642 | winner. Remember the best variable, and the direction for x<i> and x<k>. |
---|
643 | */ |
---|
644 | // try down |
---|
645 | for (k = i + 1; k < endInner; k++) { |
---|
646 | if (!maxTries) |
---|
647 | break; |
---|
648 | maxTries--; |
---|
649 | if ((way[k]&1) != 0) { |
---|
650 | // try down |
---|
651 | if (-objectiveCoefficient - cost[k] < bestChange) { |
---|
652 | // see if feasible down |
---|
653 | bool good = true; |
---|
654 | int numberMarked = 0; |
---|
655 | int kColumn = integerVariable[k]; |
---|
656 | for (j = columnStart[kColumn]; |
---|
657 | j < columnStart[kColumn] + columnLength[kColumn]; j++) { |
---|
658 | int iRow = row[j]; |
---|
659 | double newValue = rowActivity[iRow] - element[j]; |
---|
660 | if (newValue < rowLower[iRow] - primalTolerance || |
---|
661 | newValue > rowUpper[iRow] + primalTolerance) { |
---|
662 | good = false; |
---|
663 | break; |
---|
664 | } else if (mark[iRow]) { |
---|
665 | // made feasible |
---|
666 | numberMarked++; |
---|
667 | } |
---|
668 | } |
---|
669 | if (good && numberMarked == numberInfeasible) { |
---|
670 | // better solution |
---|
671 | goodK = k; |
---|
672 | wayK = -1; |
---|
673 | wayI = -1; |
---|
674 | bestChange = -objectiveCoefficient - cost[k]; |
---|
675 | } |
---|
676 | } |
---|
677 | } |
---|
678 | if ((way[k]&2) != 0) { |
---|
679 | // try up |
---|
680 | if (-objectiveCoefficient + cost[k] < bestChange) { |
---|
681 | // see if feasible up |
---|
682 | bool good = true; |
---|
683 | int numberMarked = 0; |
---|
684 | int kColumn = integerVariable[k]; |
---|
685 | for (j = columnStart[kColumn]; |
---|
686 | j < columnStart[kColumn] + columnLength[kColumn]; j++) { |
---|
687 | int iRow = row[j]; |
---|
688 | double newValue = rowActivity[iRow] + element[j]; |
---|
689 | if (newValue < rowLower[iRow] - primalTolerance || |
---|
690 | newValue > rowUpper[iRow] + primalTolerance) { |
---|
691 | good = false; |
---|
692 | break; |
---|
693 | } else if (mark[iRow]) { |
---|
694 | // made feasible |
---|
695 | numberMarked++; |
---|
696 | } |
---|
697 | } |
---|
698 | if (good && numberMarked == numberInfeasible) { |
---|
699 | // better solution |
---|
700 | goodK = k; |
---|
701 | wayK = 1; |
---|
702 | wayI = -1; |
---|
703 | bestChange = -objectiveCoefficient + cost[k]; |
---|
704 | } |
---|
705 | } |
---|
706 | } |
---|
707 | } |
---|
708 | /* |
---|
709 | Remove effect of decrementing x<i> by restoring original lhs values. |
---|
710 | */ |
---|
711 | // restore row activities |
---|
712 | for (j = columnStart[iColumn]; |
---|
713 | j < columnStart[iColumn] + columnLength[iColumn]; j++) { |
---|
714 | int iRow = row[j]; |
---|
715 | rowActivity[iRow] = save[iRow]; |
---|
716 | mark[iRow] = 0; |
---|
717 | } |
---|
718 | } |
---|
719 | /* |
---|
720 | Try to increment x<i>. Actions as for decrement. |
---|
721 | */ |
---|
722 | if ((way[i]&2) != 0) { |
---|
723 | int numberInfeasible = 0; |
---|
724 | // save row activities and adjust |
---|
725 | for (j = columnStart[iColumn]; |
---|
726 | j < columnStart[iColumn] + columnLength[iColumn]; j++) { |
---|
727 | int iRow = row[j]; |
---|
728 | save[iRow] = rowActivity[iRow]; |
---|
729 | rowActivity[iRow] += element[j]; |
---|
730 | if (rowActivity[iRow] < rowLower[iRow] - primalTolerance || |
---|
731 | rowActivity[iRow] > rowUpper[iRow] + primalTolerance) { |
---|
732 | // mark row |
---|
733 | mark[iRow] = 1; |
---|
734 | numberInfeasible++; |
---|
735 | } |
---|
736 | } |
---|
737 | // try up |
---|
738 | for (k = i + 1; k < endInner; k++) { |
---|
739 | if (!maxTries) |
---|
740 | break; |
---|
741 | if ((way[k]&1) != 0) { |
---|
742 | // try down |
---|
743 | if (objectiveCoefficient - cost[k] < bestChange) { |
---|
744 | // see if feasible down |
---|
745 | bool good = true; |
---|
746 | int numberMarked = 0; |
---|
747 | int kColumn = integerVariable[k]; |
---|
748 | for (j = columnStart[kColumn]; |
---|
749 | j < columnStart[kColumn] + columnLength[kColumn]; j++) { |
---|
750 | int iRow = row[j]; |
---|
751 | double newValue = rowActivity[iRow] - element[j]; |
---|
752 | if (newValue < rowLower[iRow] - primalTolerance || |
---|
753 | newValue > rowUpper[iRow] + primalTolerance) { |
---|
754 | good = false; |
---|
755 | break; |
---|
756 | } else if (mark[iRow]) { |
---|
757 | // made feasible |
---|
758 | numberMarked++; |
---|
759 | } |
---|
760 | } |
---|
761 | if (good && numberMarked == numberInfeasible) { |
---|
762 | // better solution |
---|
763 | goodK = k; |
---|
764 | wayK = -1; |
---|
765 | wayI = 1; |
---|
766 | bestChange = objectiveCoefficient - cost[k]; |
---|
767 | } |
---|
768 | } |
---|
769 | } |
---|
770 | if ((way[k]&2) != 0) { |
---|
771 | // try up |
---|
772 | if (objectiveCoefficient + cost[k] < bestChange) { |
---|
773 | // see if feasible up |
---|
774 | bool good = true; |
---|
775 | int numberMarked = 0; |
---|
776 | int kColumn = integerVariable[k]; |
---|
777 | for (j = columnStart[kColumn]; |
---|
778 | j < columnStart[kColumn] + columnLength[kColumn]; j++) { |
---|
779 | int iRow = row[j]; |
---|
780 | double newValue = rowActivity[iRow] + element[j]; |
---|
781 | if (newValue < rowLower[iRow] - primalTolerance || |
---|
782 | newValue > rowUpper[iRow] + primalTolerance) { |
---|
783 | good = false; |
---|
784 | break; |
---|
785 | } else if (mark[iRow]) { |
---|
786 | // made feasible |
---|
787 | numberMarked++; |
---|
788 | } |
---|
789 | } |
---|
790 | if (good && numberMarked == numberInfeasible) { |
---|
791 | // better solution |
---|
792 | goodK = k; |
---|
793 | wayK = 1; |
---|
794 | wayI = 1; |
---|
795 | bestChange = objectiveCoefficient + cost[k]; |
---|
796 | } |
---|
797 | } |
---|
798 | } |
---|
799 | } |
---|
800 | // restore row activities |
---|
801 | for (j = columnStart[iColumn]; |
---|
802 | j < columnStart[iColumn] + columnLength[iColumn]; j++) { |
---|
803 | int iRow = row[j]; |
---|
804 | rowActivity[iRow] = save[iRow]; |
---|
805 | mark[iRow] = 0; |
---|
806 | } |
---|
807 | } |
---|
808 | /* |
---|
809 | We've found a pair x<i> and x<k> which produce a better solution. Update our |
---|
810 | notion of current solution to match. |
---|
811 | |
---|
812 | Why does this not update newSolutionValue? |
---|
813 | */ |
---|
814 | if (goodK >= 0) { |
---|
815 | // we found something - update solution |
---|
816 | for (j = columnStart[iColumn]; |
---|
817 | j < columnStart[iColumn] + columnLength[iColumn]; j++) { |
---|
818 | int iRow = row[j]; |
---|
819 | rowActivity[iRow] += wayI * element[j]; |
---|
820 | } |
---|
821 | newSolution[iColumn] += wayI; |
---|
822 | int kColumn = integerVariable[goodK]; |
---|
823 | for (j = columnStart[kColumn]; |
---|
824 | j < columnStart[kColumn] + columnLength[kColumn]; j++) { |
---|
825 | int iRow = row[j]; |
---|
826 | rowActivity[iRow] += wayK * element[j]; |
---|
827 | } |
---|
828 | newSolution[kColumn] += wayK; |
---|
829 | /* |
---|
830 | Adjust motion range for x<k>. We may have banged up against a bound with that |
---|
831 | last move. |
---|
832 | */ |
---|
833 | // See if k can go further ? |
---|
834 | const OsiObject * object = model_->object(goodK); |
---|
835 | // get original bounds |
---|
836 | double originalLower; |
---|
837 | double originalUpper; |
---|
838 | getIntegerInformation( object, originalLower, originalUpper); |
---|
839 | |
---|
840 | double value = newSolution[kColumn]; |
---|
841 | int iway = 0; |
---|
842 | |
---|
843 | if (value > originalLower + 0.5) |
---|
844 | iway = 1; |
---|
845 | if (value < originalUpper - 0.5) |
---|
846 | iway |= 2; |
---|
847 | way[goodK] = static_cast<char>(iway); |
---|
848 | totalChange += bestChange; |
---|
849 | } |
---|
850 | } |
---|
851 | /* |
---|
852 | End of loop to try increment/decrement of integer variables. |
---|
853 | |
---|
854 | newSolutionValue does not necessarily match the current newSolution, and |
---|
855 | bestChange simply reflects the best single change. Still, that's sufficient |
---|
856 | to indicate that there's been at least one change. Check that we really do |
---|
857 | have a valid solution. |
---|
858 | */ |
---|
859 | if (totalChange + newSolutionValue < solutionValue) { |
---|
860 | // paranoid check |
---|
861 | memset(rowActivity, 0, numberRows*sizeof(double)); |
---|
862 | |
---|
863 | for (i = 0; i < numberColumns; i++) { |
---|
864 | int j; |
---|
865 | double value = newSolution[i]; |
---|
866 | if (value) { |
---|
867 | for (j = columnStart[i]; |
---|
868 | j < columnStart[i] + columnLength[i]; j++) { |
---|
869 | int iRow = row[j]; |
---|
870 | rowActivity[iRow] += value * element[j]; |
---|
871 | } |
---|
872 | } |
---|
873 | } |
---|
874 | int numberBad = 0; |
---|
875 | double sumBad = 0.0; |
---|
876 | // check was approximately feasible |
---|
877 | for (i = 0; i < numberRows; i++) { |
---|
878 | if (rowActivity[i] < rowLower[i]) { |
---|
879 | sumBad += rowLower[i] - rowActivity[i]; |
---|
880 | if (rowActivity[i] < rowLower[i] - 10.0*primalTolerance) |
---|
881 | numberBad++; |
---|
882 | } else if (rowActivity[i] > rowUpper[i]) { |
---|
883 | sumBad += rowUpper[i] - rowActivity[i]; |
---|
884 | if (rowActivity[i] > rowUpper[i] + 10.0*primalTolerance) |
---|
885 | numberBad++; |
---|
886 | } |
---|
887 | } |
---|
888 | if (!numberBad) { |
---|
889 | for (i = 0; i < numberIntegers; i++) { |
---|
890 | int iColumn = integerVariable[i]; |
---|
891 | const OsiObject * object = model_->object(i); |
---|
892 | // get original bounds |
---|
893 | double originalLower; |
---|
894 | double originalUpper; |
---|
895 | getIntegerInformation( object, originalLower, originalUpper); |
---|
896 | |
---|
897 | double value = newSolution[iColumn]; |
---|
898 | // if away from lower bound mark that fact |
---|
899 | if (value > originalLower) { |
---|
900 | used_[iColumn] = numberSolutions_; |
---|
901 | } |
---|
902 | } |
---|
903 | /* |
---|
904 | Copy the solution to the array returned to the client. Grab a basis from |
---|
905 | the solver (which, if it exists, is almost certainly infeasible, but it |
---|
906 | should be ok for a dual start). The value returned as solutionValue is |
---|
907 | conservative because of handling of newSolutionValue and bestChange, as |
---|
908 | described above. |
---|
909 | */ |
---|
910 | // new solution |
---|
911 | memcpy(betterSolution, newSolution, numberColumns*sizeof(double)); |
---|
912 | CoinWarmStartBasis * basis = |
---|
913 | dynamic_cast<CoinWarmStartBasis *>(solver->getWarmStart()) ; |
---|
914 | if (basis) { |
---|
915 | model_->setBestSolutionBasis(* basis); |
---|
916 | delete basis; |
---|
917 | } |
---|
918 | returnCode = 1; |
---|
919 | solutionValue = newSolutionValue + bestChange; |
---|
920 | } else { |
---|
921 | // bad solution - should not happen so debug if see message |
---|
922 | COIN_DETAIL_PRINT(printf("Local search got bad solution with %d infeasibilities summing to %g\n", |
---|
923 | numberBad, sumBad)); |
---|
924 | } |
---|
925 | } |
---|
926 | // This is just a copy! |
---|
927 | delete [] integerVariable; |
---|
928 | } |
---|
929 | /* |
---|
930 | We're done. Clean up. |
---|
931 | */ |
---|
932 | delete [] newSolution; |
---|
933 | delete [] rowActivity; |
---|
934 | delete [] way; |
---|
935 | delete [] cost; |
---|
936 | delete [] save; |
---|
937 | delete [] mark; |
---|
938 | /* |
---|
939 | Do we want to try swapping values between solutions? |
---|
940 | swap_ is set elsewhere; it's not adjusted during heuristic execution. |
---|
941 | |
---|
942 | Again, redundant test. We shouldn't be here if numberSolutions_ = 1. |
---|
943 | */ |
---|
944 | if (numberSolutions_ > 1 && (swap%10) == 1) { |
---|
945 | // try merge |
---|
946 | int returnCode2 = solutionFix( solutionValue, betterSolution, NULL); |
---|
947 | if (returnCode2) |
---|
948 | returnCode = 1; |
---|
949 | } |
---|
950 | return returnCode; |
---|
951 | } |
---|
952 | // update model |
---|
953 | void CbcHeuristicLocal::setModel(CbcModel * model) |
---|
954 | { |
---|
955 | model_ = model; |
---|
956 | // Get a copy of original matrix |
---|
957 | assert(model_->solver()); |
---|
958 | if (model_->solver()->getNumRows()) { |
---|
959 | matrix_ = *model_->solver()->getMatrixByCol(); |
---|
960 | } |
---|
961 | delete [] used_; |
---|
962 | int numberColumns = model->solver()->getNumCols(); |
---|
963 | used_ = new int[numberColumns]; |
---|
964 | memset(used_, 0, numberColumns*sizeof(int)); |
---|
965 | } |
---|
966 | |
---|
967 | // Default Constructor |
---|
968 | CbcHeuristicProximity::CbcHeuristicProximity() |
---|
969 | : CbcHeuristic() |
---|
970 | { |
---|
971 | increment_ = 0.01; |
---|
972 | feasibilityPump_ = NULL; |
---|
973 | numberSolutions_ = 0; |
---|
974 | used_ = NULL; |
---|
975 | lastRunDeep_ = -1000000; |
---|
976 | switches_ |= 16; // needs a new solution |
---|
977 | } |
---|
978 | |
---|
979 | // Constructor with model - assumed before cuts |
---|
980 | |
---|
981 | CbcHeuristicProximity::CbcHeuristicProximity(CbcModel & model) |
---|
982 | : CbcHeuristic(model) |
---|
983 | { |
---|
984 | increment_ = 0.01; |
---|
985 | feasibilityPump_ = NULL; |
---|
986 | numberSolutions_ = 0; |
---|
987 | lastRunDeep_ = -1000000; |
---|
988 | switches_ |= 16; // needs a new solution |
---|
989 | int numberColumns = model.solver()->getNumCols(); |
---|
990 | used_ = new int[numberColumns]; |
---|
991 | memset(used_, 0, numberColumns*sizeof(int)); |
---|
992 | } |
---|
993 | |
---|
994 | // Destructor |
---|
995 | CbcHeuristicProximity::~CbcHeuristicProximity () |
---|
996 | { |
---|
997 | delete feasibilityPump_; |
---|
998 | delete [] used_; |
---|
999 | } |
---|
1000 | |
---|
1001 | // Clone |
---|
1002 | CbcHeuristic * |
---|
1003 | CbcHeuristicProximity::clone() const |
---|
1004 | { |
---|
1005 | return new CbcHeuristicProximity(*this); |
---|
1006 | } |
---|
1007 | // Create C++ lines to get to current state |
---|
1008 | void |
---|
1009 | CbcHeuristicProximity::generateCpp( FILE * fp) |
---|
1010 | { |
---|
1011 | CbcHeuristicProximity other; |
---|
1012 | fprintf(fp, "0#include \"CbcHeuristicProximity.hpp\"\n"); |
---|
1013 | fprintf(fp, "3 CbcHeuristicProximity heuristicProximity(*cbcModel);\n"); |
---|
1014 | CbcHeuristic::generateCpp(fp, "heuristicProximity"); |
---|
1015 | fprintf(fp, "3 cbcModel->addHeuristic(&heuristicProximity);\n"); |
---|
1016 | } |
---|
1017 | |
---|
1018 | // Copy constructor |
---|
1019 | CbcHeuristicProximity::CbcHeuristicProximity(const CbcHeuristicProximity & rhs) |
---|
1020 | : |
---|
1021 | CbcHeuristic(rhs), |
---|
1022 | numberSolutions_(rhs.numberSolutions_) |
---|
1023 | { |
---|
1024 | increment_ = rhs.increment_; |
---|
1025 | feasibilityPump_ = NULL; |
---|
1026 | if (model_ && rhs.used_) { |
---|
1027 | int numberColumns = model_->solver()->getNumCols(); |
---|
1028 | used_ = CoinCopyOfArray(rhs.used_, numberColumns); |
---|
1029 | if (rhs.feasibilityPump_) |
---|
1030 | feasibilityPump_ = new CbcHeuristicFPump(*rhs.feasibilityPump_); |
---|
1031 | } else { |
---|
1032 | used_ = NULL; |
---|
1033 | } |
---|
1034 | } |
---|
1035 | |
---|
1036 | // Assignment operator |
---|
1037 | CbcHeuristicProximity & |
---|
1038 | CbcHeuristicProximity::operator=( const CbcHeuristicProximity & rhs) |
---|
1039 | { |
---|
1040 | if (this != &rhs) { |
---|
1041 | CbcHeuristic::operator=(rhs); |
---|
1042 | increment_ = rhs.increment_; |
---|
1043 | numberSolutions_ = rhs.numberSolutions_; |
---|
1044 | delete [] used_; |
---|
1045 | delete feasibilityPump_; |
---|
1046 | feasibilityPump_ = NULL; |
---|
1047 | if (model_ && rhs.used_) { |
---|
1048 | int numberColumns = model_->solver()->getNumCols(); |
---|
1049 | used_ = CoinCopyOfArray(rhs.used_, numberColumns); |
---|
1050 | if (rhs.feasibilityPump_) |
---|
1051 | feasibilityPump_ = new CbcHeuristicFPump(*rhs.feasibilityPump_); |
---|
1052 | } else { |
---|
1053 | used_ = NULL; |
---|
1054 | } |
---|
1055 | } |
---|
1056 | return *this; |
---|
1057 | } |
---|
1058 | |
---|
1059 | // Resets stuff if model changes |
---|
1060 | void |
---|
1061 | CbcHeuristicProximity::resetModel(CbcModel * /*model*/) |
---|
1062 | { |
---|
1063 | //CbcHeuristic::resetModel(model); |
---|
1064 | delete [] used_; |
---|
1065 | if (model_ && used_) { |
---|
1066 | int numberColumns = model_->solver()->getNumCols(); |
---|
1067 | used_ = new int[numberColumns]; |
---|
1068 | memset(used_, 0, numberColumns*sizeof(int)); |
---|
1069 | } else { |
---|
1070 | used_ = NULL; |
---|
1071 | } |
---|
1072 | } |
---|
1073 | /* |
---|
1074 | Run a mini-BaB search after changing objective |
---|
1075 | |
---|
1076 | Return values are: |
---|
1077 | 1: smallBranchAndBound found a solution |
---|
1078 | 0: everything else |
---|
1079 | |
---|
1080 | The degree of overload as return codes from smallBranchAndBound are folded |
---|
1081 | into 0 is such that it's impossible to distinguish return codes that really |
---|
1082 | require attention from a simple `nothing of interest'. |
---|
1083 | */ |
---|
1084 | int |
---|
1085 | CbcHeuristicProximity::solution(double & solutionValue, |
---|
1086 | double * betterSolution) |
---|
1087 | { |
---|
1088 | if (feasibilityPumpOptions_ == -3 && numCouldRun_==0 && |
---|
1089 | !feasibilityPump_ ) { |
---|
1090 | // clone feasibility pump |
---|
1091 | for (int i = 0; i < model_->numberHeuristics(); i++) { |
---|
1092 | const CbcHeuristicFPump* pump = |
---|
1093 | dynamic_cast<const CbcHeuristicFPump*>(model_->heuristic(i)); |
---|
1094 | if (pump) { |
---|
1095 | feasibilityPump_ = new CbcHeuristicFPump(*pump); |
---|
1096 | break; |
---|
1097 | } |
---|
1098 | } |
---|
1099 | } |
---|
1100 | /* |
---|
1101 | Execute only if a new solution has been discovered since the last time we |
---|
1102 | were called. |
---|
1103 | */ |
---|
1104 | |
---|
1105 | numCouldRun_++; |
---|
1106 | int nodeCount = model_->getNodeCount(); |
---|
1107 | if (numberSolutions_ == model_->getSolutionCount()) |
---|
1108 | return 0; |
---|
1109 | if (!model_->bestSolution()) |
---|
1110 | return 0; // odd - because in parallel mode |
---|
1111 | numberSolutions_ = model_->getSolutionCount(); |
---|
1112 | lastRunDeep_ = nodeCount; |
---|
1113 | numRuns_++; |
---|
1114 | //howOftenShallow_ = numberSolutions_; |
---|
1115 | |
---|
1116 | /* |
---|
1117 | Load up a new solver with the solution. |
---|
1118 | |
---|
1119 | Why continuousSolver(), as opposed to solver()? |
---|
1120 | */ |
---|
1121 | OsiSolverInterface * newSolver = model_->continuousSolver()->clone(); |
---|
1122 | int numberColumns=newSolver->getNumCols(); |
---|
1123 | double * obj = CoinCopyOfArray(newSolver->getObjCoefficients(),numberColumns); |
---|
1124 | int * indices = new int [numberColumns]; |
---|
1125 | int n=0; |
---|
1126 | for (int i=0;i<numberColumns;i++) { |
---|
1127 | if (obj[i]) { |
---|
1128 | indices[n]=i; |
---|
1129 | obj[n++]=obj[i]; |
---|
1130 | } |
---|
1131 | } |
---|
1132 | double cutoff=model_->getCutoff(); |
---|
1133 | assert (cutoff<1.0e20); |
---|
1134 | if (model_->getCutoffIncrement()<1.0e-4) { |
---|
1135 | cutoff -= increment_; |
---|
1136 | } |
---|
1137 | double offset; |
---|
1138 | newSolver->getDblParam(OsiObjOffset, offset); |
---|
1139 | newSolver->setDblParam(OsiObjOffset, 0.0); |
---|
1140 | newSolver->addRow(n,indices,obj,-COIN_DBL_MAX,cutoff+offset); |
---|
1141 | delete [] indices; |
---|
1142 | memset(obj,0,numberColumns*sizeof(double)); |
---|
1143 | newSolver->setDblParam(OsiDualObjectiveLimit, 1.0e20); |
---|
1144 | int numberIntegers = model_->numberIntegers(); |
---|
1145 | const int * integerVariable = model_->integerVariable(); |
---|
1146 | const double * solutionIn = model_->bestSolution(); |
---|
1147 | for (int i = 0; i < numberIntegers; i++) { |
---|
1148 | int iColumn = integerVariable[i]; |
---|
1149 | if (fabs(solutionIn[iColumn])<1.0e-5) |
---|
1150 | obj[iColumn]=1.0; |
---|
1151 | else if (fabs(solutionIn[iColumn]-1.0)<1.0e-5) |
---|
1152 | obj[iColumn]=-1.0; |
---|
1153 | } |
---|
1154 | newSolver->setObjective(obj); |
---|
1155 | delete [] obj; |
---|
1156 | //newSolver->writeMps("xxxx"); |
---|
1157 | int maxSolutions = model_->getMaximumSolutions(); |
---|
1158 | model_->setMaximumSolutions(1); |
---|
1159 | bool pumpAdded = false; |
---|
1160 | if (feasibilityPumpOptions_ == -3 && feasibilityPump_) { |
---|
1161 | // add back feasibility pump |
---|
1162 | pumpAdded = true; |
---|
1163 | for (int i = 0; i < model_->numberHeuristics(); i++) { |
---|
1164 | const CbcHeuristicFPump* pump = |
---|
1165 | dynamic_cast<const CbcHeuristicFPump*>(model_->heuristic(i)); |
---|
1166 | if (pump) { |
---|
1167 | pumpAdded = false; |
---|
1168 | break; |
---|
1169 | } |
---|
1170 | } |
---|
1171 | if (pumpAdded) |
---|
1172 | model_->addHeuristic(feasibilityPump_); |
---|
1173 | } |
---|
1174 | int returnCode = |
---|
1175 | smallBranchAndBound(newSolver, numberNodes_, betterSolution, solutionValue, |
---|
1176 | 1.0e20, "CbcHeuristicProximity"); |
---|
1177 | if (pumpAdded) { |
---|
1178 | // take off feasibility pump |
---|
1179 | int lastHeuristic = model_->numberHeuristics()-1; |
---|
1180 | model_->setNumberHeuristics(lastHeuristic); |
---|
1181 | delete model_->heuristic(lastHeuristic); |
---|
1182 | } |
---|
1183 | model_->setMaximumSolutions(maxSolutions); |
---|
1184 | /* |
---|
1185 | -2 is return due to user event, and -1 is overloaded with what look to be |
---|
1186 | two contradictory meanings. |
---|
1187 | */ |
---|
1188 | if (returnCode < 0) { |
---|
1189 | returnCode = 0; |
---|
1190 | } |
---|
1191 | /* |
---|
1192 | If the result is complete exploration with a solution (3) or proven |
---|
1193 | infeasibility (2), we could generate a cut (the AI folks would call it a |
---|
1194 | nogood) to prevent us from going down this route in the future. |
---|
1195 | */ |
---|
1196 | if ((returnCode&2) != 0) { |
---|
1197 | // could add cut |
---|
1198 | returnCode &= ~2; |
---|
1199 | } |
---|
1200 | char proxPrint[200]; |
---|
1201 | if ((returnCode&1) != 0) { |
---|
1202 | // redo objective |
---|
1203 | const double * obj = model_->continuousSolver()->getObjCoefficients(); |
---|
1204 | solutionValue = - offset; |
---|
1205 | int sumIncrease=0.0; |
---|
1206 | int sumDecrease=0.0; |
---|
1207 | int numberIncrease=0; |
---|
1208 | int numberDecrease=0; |
---|
1209 | for (int i=0;i<numberColumns;i++) { |
---|
1210 | solutionValue += obj[i]*betterSolution[i]; |
---|
1211 | if (model_->isInteger(i)) { |
---|
1212 | int change=static_cast<int>(floor(solutionIn[i]-betterSolution[i]+0.5)); |
---|
1213 | if (change>0) { |
---|
1214 | numberIncrease++; |
---|
1215 | sumIncrease+=change; |
---|
1216 | } else if (change<0) { |
---|
1217 | numberDecrease++; |
---|
1218 | sumDecrease-=change; |
---|
1219 | } |
---|
1220 | } |
---|
1221 | } |
---|
1222 | sprintf(proxPrint,"Proximity search ran %d nodes (out of %d) - in new solution %d increased (%d), %d decreased (%d)", |
---|
1223 | numberNodesDone_,numberNodes_, |
---|
1224 | numberIncrease,sumIncrease,numberDecrease,sumDecrease); |
---|
1225 | if (!numberIncrease&&!numberDecrease) { |
---|
1226 | // somehow tolerances are such that we can slip through |
---|
1227 | // change for next time |
---|
1228 | increment_ += CoinMax(increment_,fabs(solutionValue+offset)*1.0e-10); |
---|
1229 | } |
---|
1230 | } else { |
---|
1231 | sprintf(proxPrint,"Proximity search ran %d nodes - no new solution", |
---|
1232 | numberNodesDone_); |
---|
1233 | } |
---|
1234 | model_->messageHandler()->message(CBC_FPUMP1, model_->messages()) |
---|
1235 | << proxPrint |
---|
1236 | << CoinMessageEol; |
---|
1237 | |
---|
1238 | delete newSolver; |
---|
1239 | return returnCode; |
---|
1240 | } |
---|
1241 | // update model |
---|
1242 | void CbcHeuristicProximity::setModel(CbcModel * model) |
---|
1243 | { |
---|
1244 | model_ = model; |
---|
1245 | // Get a copy of original matrix |
---|
1246 | assert(model_->solver()); |
---|
1247 | delete [] used_; |
---|
1248 | int numberColumns = model->solver()->getNumCols(); |
---|
1249 | used_ = new int[numberColumns]; |
---|
1250 | memset(used_, 0, numberColumns*sizeof(int)); |
---|
1251 | } |
---|
1252 | |
---|
1253 | // Default Constructor |
---|
1254 | CbcHeuristicNaive::CbcHeuristicNaive() |
---|
1255 | : CbcHeuristic() |
---|
1256 | { |
---|
1257 | large_ = 1.0e6; |
---|
1258 | } |
---|
1259 | |
---|
1260 | // Constructor with model - assumed before cuts |
---|
1261 | |
---|
1262 | CbcHeuristicNaive::CbcHeuristicNaive(CbcModel & model) |
---|
1263 | : CbcHeuristic(model) |
---|
1264 | { |
---|
1265 | large_ = 1.0e6; |
---|
1266 | } |
---|
1267 | |
---|
1268 | // Destructor |
---|
1269 | CbcHeuristicNaive::~CbcHeuristicNaive () |
---|
1270 | { |
---|
1271 | } |
---|
1272 | |
---|
1273 | // Clone |
---|
1274 | CbcHeuristic * |
---|
1275 | CbcHeuristicNaive::clone() const |
---|
1276 | { |
---|
1277 | return new CbcHeuristicNaive(*this); |
---|
1278 | } |
---|
1279 | // Create C++ lines to get to current state |
---|
1280 | void |
---|
1281 | CbcHeuristicNaive::generateCpp( FILE * fp) |
---|
1282 | { |
---|
1283 | CbcHeuristicNaive other; |
---|
1284 | fprintf(fp, "0#include \"CbcHeuristicProximity.hpp\"\n"); |
---|
1285 | fprintf(fp, "3 CbcHeuristicNaive naive(*cbcModel);\n"); |
---|
1286 | CbcHeuristic::generateCpp(fp, "naive"); |
---|
1287 | if (large_ != other.large_) |
---|
1288 | fprintf(fp, "3 naive.setLarge(%g);\n", large_); |
---|
1289 | else |
---|
1290 | fprintf(fp, "4 naive.setLarge(%g);\n", large_); |
---|
1291 | fprintf(fp, "3 cbcModel->addHeuristic(&naive);\n"); |
---|
1292 | } |
---|
1293 | |
---|
1294 | // Copy constructor |
---|
1295 | CbcHeuristicNaive::CbcHeuristicNaive(const CbcHeuristicNaive & rhs) |
---|
1296 | : |
---|
1297 | CbcHeuristic(rhs), |
---|
1298 | large_(rhs.large_) |
---|
1299 | { |
---|
1300 | } |
---|
1301 | |
---|
1302 | // Assignment operator |
---|
1303 | CbcHeuristicNaive & |
---|
1304 | CbcHeuristicNaive::operator=( const CbcHeuristicNaive & rhs) |
---|
1305 | { |
---|
1306 | if (this != &rhs) { |
---|
1307 | CbcHeuristic::operator=(rhs); |
---|
1308 | large_ = rhs.large_; |
---|
1309 | } |
---|
1310 | return *this; |
---|
1311 | } |
---|
1312 | |
---|
1313 | // Resets stuff if model changes |
---|
1314 | void |
---|
1315 | CbcHeuristicNaive::resetModel(CbcModel * model) |
---|
1316 | { |
---|
1317 | CbcHeuristic::resetModel(model); |
---|
1318 | } |
---|
1319 | int |
---|
1320 | CbcHeuristicNaive::solution(double & solutionValue, |
---|
1321 | double * betterSolution) |
---|
1322 | { |
---|
1323 | numCouldRun_++; |
---|
1324 | // See if to do |
---|
1325 | bool atRoot = model_->getNodeCount() == 0; |
---|
1326 | int passNumber = model_->getCurrentPassNumber(); |
---|
1327 | if (!when() || (when() == 1 && model_->phase() != 1) || !atRoot || passNumber > 1) |
---|
1328 | return 0; // switched off |
---|
1329 | // Don't do if it was this heuristic which found solution! |
---|
1330 | if (this == model_->lastHeuristic()) |
---|
1331 | return 0; |
---|
1332 | numRuns_++; |
---|
1333 | double cutoff; |
---|
1334 | model_->solver()->getDblParam(OsiDualObjectiveLimit, cutoff); |
---|
1335 | double direction = model_->solver()->getObjSense(); |
---|
1336 | cutoff *= direction; |
---|
1337 | cutoff = CoinMin(cutoff, solutionValue); |
---|
1338 | OsiSolverInterface * solver = model_->continuousSolver(); |
---|
1339 | if (!solver) |
---|
1340 | solver = model_->solver(); |
---|
1341 | const double * colLower = solver->getColLower(); |
---|
1342 | const double * colUpper = solver->getColUpper(); |
---|
1343 | const double * objective = solver->getObjCoefficients(); |
---|
1344 | |
---|
1345 | int numberColumns = model_->getNumCols(); |
---|
1346 | int numberIntegers = model_->numberIntegers(); |
---|
1347 | const int * integerVariable = model_->integerVariable(); |
---|
1348 | |
---|
1349 | int i; |
---|
1350 | bool solutionFound = false; |
---|
1351 | CoinWarmStartBasis saveBasis; |
---|
1352 | CoinWarmStartBasis * basis = |
---|
1353 | dynamic_cast<CoinWarmStartBasis *>(solver->getWarmStart()) ; |
---|
1354 | if (basis) { |
---|
1355 | saveBasis = * basis; |
---|
1356 | delete basis; |
---|
1357 | } |
---|
1358 | // First just fix all integers as close to zero as possible |
---|
1359 | OsiSolverInterface * newSolver = cloneBut(7); // wassolver->clone(); |
---|
1360 | for (i = 0; i < numberIntegers; i++) { |
---|
1361 | int iColumn = integerVariable[i]; |
---|
1362 | double lower = colLower[iColumn]; |
---|
1363 | double upper = colUpper[iColumn]; |
---|
1364 | double value; |
---|
1365 | if (lower > 0.0) |
---|
1366 | value = lower; |
---|
1367 | else if (upper < 0.0) |
---|
1368 | value = upper; |
---|
1369 | else |
---|
1370 | value = 0.0; |
---|
1371 | newSolver->setColLower(iColumn, value); |
---|
1372 | newSolver->setColUpper(iColumn, value); |
---|
1373 | } |
---|
1374 | newSolver->initialSolve(); |
---|
1375 | if (newSolver->isProvenOptimal()) { |
---|
1376 | double solValue = newSolver->getObjValue() * direction ; |
---|
1377 | if (solValue < cutoff) { |
---|
1378 | // we have a solution |
---|
1379 | solutionFound = true; |
---|
1380 | solutionValue = solValue; |
---|
1381 | memcpy(betterSolution, newSolver->getColSolution(), |
---|
1382 | numberColumns*sizeof(double)); |
---|
1383 | COIN_DETAIL_PRINT(printf("Naive fixing close to zero gave solution of %g\n", solutionValue)); |
---|
1384 | cutoff = solValue - model_->getCutoffIncrement(); |
---|
1385 | } |
---|
1386 | } |
---|
1387 | // Now fix all integers as close to zero if not zero or large cost |
---|
1388 | int nFix = 0; |
---|
1389 | for (i = 0; i < numberIntegers; i++) { |
---|
1390 | int iColumn = integerVariable[i]; |
---|
1391 | double lower = colLower[iColumn]; |
---|
1392 | double upper = colUpper[iColumn]; |
---|
1393 | double value; |
---|
1394 | if (fabs(objective[i]) > 0.0 && fabs(objective[i]) < large_) { |
---|
1395 | nFix++; |
---|
1396 | if (lower > 0.0) |
---|
1397 | value = lower; |
---|
1398 | else if (upper < 0.0) |
---|
1399 | value = upper; |
---|
1400 | else |
---|
1401 | value = 0.0; |
---|
1402 | newSolver->setColLower(iColumn, value); |
---|
1403 | newSolver->setColUpper(iColumn, value); |
---|
1404 | } else { |
---|
1405 | // set back to original |
---|
1406 | newSolver->setColLower(iColumn, lower); |
---|
1407 | newSolver->setColUpper(iColumn, upper); |
---|
1408 | } |
---|
1409 | } |
---|
1410 | const double * solution = solver->getColSolution(); |
---|
1411 | if (nFix) { |
---|
1412 | newSolver->setWarmStart(&saveBasis); |
---|
1413 | newSolver->setColSolution(solution); |
---|
1414 | newSolver->initialSolve(); |
---|
1415 | if (newSolver->isProvenOptimal()) { |
---|
1416 | double solValue = newSolver->getObjValue() * direction ; |
---|
1417 | if (solValue < cutoff) { |
---|
1418 | // try branch and bound |
---|
1419 | double * newSolution = new double [numberColumns]; |
---|
1420 | COIN_DETAIL_PRINT(printf("%d fixed after fixing costs\n", nFix)); |
---|
1421 | int returnCode = smallBranchAndBound(newSolver, |
---|
1422 | numberNodes_, newSolution, |
---|
1423 | solutionValue, |
---|
1424 | solutionValue, "CbcHeuristicNaive1"); |
---|
1425 | if (returnCode < 0) |
---|
1426 | returnCode = 0; // returned on size |
---|
1427 | if ((returnCode&2) != 0) { |
---|
1428 | // could add cut |
---|
1429 | returnCode &= ~2; |
---|
1430 | } |
---|
1431 | if (returnCode == 1) { |
---|
1432 | // solution |
---|
1433 | solutionFound = true; |
---|
1434 | memcpy(betterSolution, newSolution, |
---|
1435 | numberColumns*sizeof(double)); |
---|
1436 | COIN_DETAIL_PRINT(printf("Naive fixing zeros gave solution of %g\n", solutionValue)); |
---|
1437 | cutoff = solutionValue - model_->getCutoffIncrement(); |
---|
1438 | } |
---|
1439 | delete [] newSolution; |
---|
1440 | } |
---|
1441 | } |
---|
1442 | } |
---|
1443 | #if 1 |
---|
1444 | newSolver->setObjSense(-direction); // maximize |
---|
1445 | newSolver->setWarmStart(&saveBasis); |
---|
1446 | newSolver->setColSolution(solution); |
---|
1447 | for (int iColumn = 0; iColumn < numberColumns; iColumn++) { |
---|
1448 | double value = solution[iColumn]; |
---|
1449 | double lower = colLower[iColumn]; |
---|
1450 | double upper = colUpper[iColumn]; |
---|
1451 | double newLower; |
---|
1452 | double newUpper; |
---|
1453 | if (newSolver->isInteger(iColumn)) { |
---|
1454 | newLower = CoinMax(lower, floor(value) - 2.0); |
---|
1455 | newUpper = CoinMin(upper, ceil(value) + 2.0); |
---|
1456 | } else { |
---|
1457 | newLower = CoinMax(lower, value - 1.0e5); |
---|
1458 | newUpper = CoinMin(upper, value + 1.0e-5); |
---|
1459 | } |
---|
1460 | newSolver->setColLower(iColumn, newLower); |
---|
1461 | newSolver->setColUpper(iColumn, newUpper); |
---|
1462 | } |
---|
1463 | newSolver->initialSolve(); |
---|
1464 | if (newSolver->isProvenOptimal()) { |
---|
1465 | double solValue = newSolver->getObjValue() * direction ; |
---|
1466 | if (solValue < cutoff) { |
---|
1467 | nFix = 0; |
---|
1468 | newSolver->setObjSense(direction); // correct direction |
---|
1469 | //const double * thisSolution = newSolver->getColSolution(); |
---|
1470 | for (int iColumn = 0; iColumn < numberColumns; iColumn++) { |
---|
1471 | double value = solution[iColumn]; |
---|
1472 | double lower = colLower[iColumn]; |
---|
1473 | double upper = colUpper[iColumn]; |
---|
1474 | double newLower = lower; |
---|
1475 | double newUpper = upper; |
---|
1476 | if (newSolver->isInteger(iColumn)) { |
---|
1477 | if (value < lower + 1.0e-6) { |
---|
1478 | nFix++; |
---|
1479 | newUpper = lower; |
---|
1480 | } else if (value > upper - 1.0e-6) { |
---|
1481 | nFix++; |
---|
1482 | newLower = upper; |
---|
1483 | } else { |
---|
1484 | newLower = CoinMax(lower, floor(value) - 2.0); |
---|
1485 | newUpper = CoinMin(upper, ceil(value) + 2.0); |
---|
1486 | } |
---|
1487 | } |
---|
1488 | newSolver->setColLower(iColumn, newLower); |
---|
1489 | newSolver->setColUpper(iColumn, newUpper); |
---|
1490 | } |
---|
1491 | // try branch and bound |
---|
1492 | double * newSolution = new double [numberColumns]; |
---|
1493 | COIN_DETAIL_PRINT(printf("%d fixed after maximizing\n", nFix)); |
---|
1494 | int returnCode = smallBranchAndBound(newSolver, |
---|
1495 | numberNodes_, newSolution, |
---|
1496 | solutionValue, |
---|
1497 | solutionValue, "CbcHeuristicNaive1"); |
---|
1498 | if (returnCode < 0) |
---|
1499 | returnCode = 0; // returned on size |
---|
1500 | if ((returnCode&2) != 0) { |
---|
1501 | // could add cut |
---|
1502 | returnCode &= ~2; |
---|
1503 | } |
---|
1504 | if (returnCode == 1) { |
---|
1505 | // solution |
---|
1506 | solutionFound = true; |
---|
1507 | memcpy(betterSolution, newSolution, |
---|
1508 | numberColumns*sizeof(double)); |
---|
1509 | COIN_DETAIL_PRINT(printf("Naive maximizing gave solution of %g\n", solutionValue)); |
---|
1510 | cutoff = solutionValue - model_->getCutoffIncrement(); |
---|
1511 | } |
---|
1512 | delete [] newSolution; |
---|
1513 | } |
---|
1514 | } |
---|
1515 | #endif |
---|
1516 | delete newSolver; |
---|
1517 | return solutionFound ? 1 : 0; |
---|
1518 | } |
---|
1519 | // update model |
---|
1520 | void CbcHeuristicNaive::setModel(CbcModel * model) |
---|
1521 | { |
---|
1522 | model_ = model; |
---|
1523 | } |
---|
1524 | // Default Constructor |
---|
1525 | CbcHeuristicCrossover::CbcHeuristicCrossover() |
---|
1526 | : CbcHeuristic(), |
---|
1527 | numberSolutions_(0), |
---|
1528 | useNumber_(3) |
---|
1529 | { |
---|
1530 | setWhen(1); |
---|
1531 | } |
---|
1532 | |
---|
1533 | // Constructor with model - assumed before cuts |
---|
1534 | |
---|
1535 | CbcHeuristicCrossover::CbcHeuristicCrossover(CbcModel & model) |
---|
1536 | : CbcHeuristic(model), |
---|
1537 | numberSolutions_(0), |
---|
1538 | useNumber_(3) |
---|
1539 | { |
---|
1540 | setWhen(1); |
---|
1541 | for (int i = 0; i < 10; i++) |
---|
1542 | random_[i] = model.randomNumberGenerator()->randomDouble(); |
---|
1543 | } |
---|
1544 | |
---|
1545 | // Destructor |
---|
1546 | CbcHeuristicCrossover::~CbcHeuristicCrossover () |
---|
1547 | { |
---|
1548 | } |
---|
1549 | |
---|
1550 | // Clone |
---|
1551 | CbcHeuristic * |
---|
1552 | CbcHeuristicCrossover::clone() const |
---|
1553 | { |
---|
1554 | return new CbcHeuristicCrossover(*this); |
---|
1555 | } |
---|
1556 | // Create C++ lines to get to current state |
---|
1557 | void |
---|
1558 | CbcHeuristicCrossover::generateCpp( FILE * fp) |
---|
1559 | { |
---|
1560 | CbcHeuristicCrossover other; |
---|
1561 | fprintf(fp, "0#include \"CbcHeuristicProximity.hpp\"\n"); |
---|
1562 | fprintf(fp, "3 CbcHeuristicCrossover crossover(*cbcModel);\n"); |
---|
1563 | CbcHeuristic::generateCpp(fp, "crossover"); |
---|
1564 | if (useNumber_ != other.useNumber_) |
---|
1565 | fprintf(fp, "3 crossover.setNumberSolutions(%d);\n", useNumber_); |
---|
1566 | else |
---|
1567 | fprintf(fp, "4 crossover.setNumberSolutions(%d);\n", useNumber_); |
---|
1568 | fprintf(fp, "3 cbcModel->addHeuristic(&crossover);\n"); |
---|
1569 | } |
---|
1570 | |
---|
1571 | // Copy constructor |
---|
1572 | CbcHeuristicCrossover::CbcHeuristicCrossover(const CbcHeuristicCrossover & rhs) |
---|
1573 | : |
---|
1574 | CbcHeuristic(rhs), |
---|
1575 | attempts_(rhs.attempts_), |
---|
1576 | numberSolutions_(rhs.numberSolutions_), |
---|
1577 | useNumber_(rhs.useNumber_) |
---|
1578 | { |
---|
1579 | memcpy(random_, rhs.random_, 10*sizeof(double)); |
---|
1580 | } |
---|
1581 | |
---|
1582 | // Assignment operator |
---|
1583 | CbcHeuristicCrossover & |
---|
1584 | CbcHeuristicCrossover::operator=( const CbcHeuristicCrossover & rhs) |
---|
1585 | { |
---|
1586 | if (this != &rhs) { |
---|
1587 | CbcHeuristic::operator=(rhs); |
---|
1588 | useNumber_ = rhs.useNumber_; |
---|
1589 | attempts_ = rhs.attempts_; |
---|
1590 | numberSolutions_ = rhs.numberSolutions_; |
---|
1591 | memcpy(random_, rhs.random_, 10*sizeof(double)); |
---|
1592 | } |
---|
1593 | return *this; |
---|
1594 | } |
---|
1595 | |
---|
1596 | // Resets stuff if model changes |
---|
1597 | void |
---|
1598 | CbcHeuristicCrossover::resetModel(CbcModel * model) |
---|
1599 | { |
---|
1600 | CbcHeuristic::resetModel(model); |
---|
1601 | } |
---|
1602 | int |
---|
1603 | CbcHeuristicCrossover::solution(double & solutionValue, |
---|
1604 | double * betterSolution) |
---|
1605 | { |
---|
1606 | if (when_ == 0) |
---|
1607 | return 0; |
---|
1608 | numCouldRun_++; |
---|
1609 | bool useBest = (numberSolutions_ != model_->getSolutionCount()); |
---|
1610 | if (!useBest && (when_ % 10) == 1) |
---|
1611 | return 0; |
---|
1612 | numberSolutions_ = model_->getSolutionCount(); |
---|
1613 | OsiSolverInterface * continuousSolver = model_->continuousSolver(); |
---|
1614 | int useNumber = CoinMin(model_->numberSavedSolutions(), useNumber_); |
---|
1615 | if (useNumber < 2 || !continuousSolver) |
---|
1616 | return 0; |
---|
1617 | // Fix later |
---|
1618 | if (!useBest) |
---|
1619 | abort(); |
---|
1620 | numRuns_++; |
---|
1621 | double cutoff; |
---|
1622 | model_->solver()->getDblParam(OsiDualObjectiveLimit, cutoff); |
---|
1623 | double direction = model_->solver()->getObjSense(); |
---|
1624 | cutoff *= direction; |
---|
1625 | cutoff = CoinMin(cutoff, solutionValue); |
---|
1626 | OsiSolverInterface * solver = cloneBut(2); |
---|
1627 | // But reset bounds |
---|
1628 | solver->setColLower(continuousSolver->getColLower()); |
---|
1629 | solver->setColUpper(continuousSolver->getColUpper()); |
---|
1630 | int numberColumns = solver->getNumCols(); |
---|
1631 | // Fixed |
---|
1632 | double * fixed = new double [numberColumns]; |
---|
1633 | for (int i = 0; i < numberColumns; i++) |
---|
1634 | fixed[i] = -COIN_DBL_MAX; |
---|
1635 | int whichSolution[10]; |
---|
1636 | for (int i = 0; i < useNumber; i++) |
---|
1637 | whichSolution[i] = i; |
---|
1638 | for (int i = 0; i < useNumber; i++) { |
---|
1639 | int k = whichSolution[i]; |
---|
1640 | const double * solution = model_->savedSolution(k); |
---|
1641 | for (int j = 0; j < numberColumns; j++) { |
---|
1642 | if (solver->isInteger(j)) { |
---|
1643 | if (fixed[j] == -COIN_DBL_MAX) |
---|
1644 | fixed[j] = floor(solution[j] + 0.5); |
---|
1645 | else if (fabs(fixed[j] - solution[j]) > 1.0e-7) |
---|
1646 | fixed[j] = COIN_DBL_MAX; |
---|
1647 | } |
---|
1648 | } |
---|
1649 | } |
---|
1650 | const double * colLower = solver->getColLower(); |
---|
1651 | for (int i = 0; i < numberColumns; i++) { |
---|
1652 | if (solver->isInteger(i)) { |
---|
1653 | double value = fixed[i]; |
---|
1654 | if (value != COIN_DBL_MAX) { |
---|
1655 | if (when_ < 10) { |
---|
1656 | solver->setColLower(i, value); |
---|
1657 | solver->setColUpper(i, value); |
---|
1658 | } else if (value == colLower[i]) { |
---|
1659 | solver->setColUpper(i, value); |
---|
1660 | } |
---|
1661 | } |
---|
1662 | } |
---|
1663 | } |
---|
1664 | int returnCode = smallBranchAndBound(solver, numberNodes_, betterSolution, |
---|
1665 | solutionValue, |
---|
1666 | solutionValue, "CbcHeuristicCrossover"); |
---|
1667 | if (returnCode < 0) |
---|
1668 | returnCode = 0; // returned on size |
---|
1669 | if ((returnCode&2) != 0) { |
---|
1670 | // could add cut |
---|
1671 | returnCode &= ~2; |
---|
1672 | } |
---|
1673 | |
---|
1674 | delete solver; |
---|
1675 | return returnCode; |
---|
1676 | } |
---|
1677 | // update model |
---|
1678 | void CbcHeuristicCrossover::setModel(CbcModel * model) |
---|
1679 | { |
---|
1680 | model_ = model; |
---|
1681 | if (model) { |
---|
1682 | for (int i = 0; i < 10; i++) |
---|
1683 | random_[i] = model->randomNumberGenerator()->randomDouble(); |
---|
1684 | } |
---|
1685 | } |
---|
1686 | |
---|
1687 | |
---|