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