1 | /* $Id: CbcHeuristicGreedy.cpp 2216 2015-10-06 16:09:32Z forrest $ */ |
---|
2 | // Copyright (C) 2005, 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 "CbcStrategy.hpp" |
---|
18 | #include "CbcHeuristicGreedy.hpp" |
---|
19 | #include "CoinSort.hpp" |
---|
20 | #include "CglPreProcess.hpp" |
---|
21 | // Default Constructor |
---|
22 | CbcHeuristicGreedyCover::CbcHeuristicGreedyCover() |
---|
23 | : CbcHeuristic() |
---|
24 | { |
---|
25 | // matrix will automatically be empty |
---|
26 | originalNumberRows_ = 0; |
---|
27 | algorithm_ = 0; |
---|
28 | numberTimes_ = 100; |
---|
29 | } |
---|
30 | |
---|
31 | // Constructor from model |
---|
32 | CbcHeuristicGreedyCover::CbcHeuristicGreedyCover(CbcModel & model) |
---|
33 | : CbcHeuristic(model) |
---|
34 | { |
---|
35 | gutsOfConstructor(&model); |
---|
36 | algorithm_ = 0; |
---|
37 | numberTimes_ = 100; |
---|
38 | whereFrom_ = 1; |
---|
39 | } |
---|
40 | |
---|
41 | // Destructor |
---|
42 | CbcHeuristicGreedyCover::~CbcHeuristicGreedyCover () |
---|
43 | { |
---|
44 | } |
---|
45 | |
---|
46 | // Clone |
---|
47 | CbcHeuristic * |
---|
48 | CbcHeuristicGreedyCover::clone() const |
---|
49 | { |
---|
50 | return new CbcHeuristicGreedyCover(*this); |
---|
51 | } |
---|
52 | // Guts of constructor from a CbcModel |
---|
53 | void |
---|
54 | CbcHeuristicGreedyCover::gutsOfConstructor(CbcModel * model) |
---|
55 | { |
---|
56 | model_ = model; |
---|
57 | // Get a copy of original matrix |
---|
58 | assert(model->solver()); |
---|
59 | if (model->solver()->getNumRows()) { |
---|
60 | matrix_ = *model->solver()->getMatrixByCol(); |
---|
61 | } |
---|
62 | originalNumberRows_ = model->solver()->getNumRows(); |
---|
63 | } |
---|
64 | // Create C++ lines to get to current state |
---|
65 | void |
---|
66 | CbcHeuristicGreedyCover::generateCpp( FILE * fp) |
---|
67 | { |
---|
68 | CbcHeuristicGreedyCover other; |
---|
69 | fprintf(fp, "0#include \"CbcHeuristicGreedy.hpp\"\n"); |
---|
70 | fprintf(fp, "3 CbcHeuristicGreedyCover heuristicGreedyCover(*cbcModel);\n"); |
---|
71 | CbcHeuristic::generateCpp(fp, "heuristicGreedyCover"); |
---|
72 | if (algorithm_ != other.algorithm_) |
---|
73 | fprintf(fp, "3 heuristicGreedyCover.setAlgorithm(%d);\n", algorithm_); |
---|
74 | else |
---|
75 | fprintf(fp, "4 heuristicGreedyCover.setAlgorithm(%d);\n", algorithm_); |
---|
76 | if (numberTimes_ != other.numberTimes_) |
---|
77 | fprintf(fp, "3 heuristicGreedyCover.setNumberTimes(%d);\n", numberTimes_); |
---|
78 | else |
---|
79 | fprintf(fp, "4 heuristicGreedyCover.setNumberTimes(%d);\n", numberTimes_); |
---|
80 | fprintf(fp, "3 cbcModel->addHeuristic(&heuristicGreedyCover);\n"); |
---|
81 | } |
---|
82 | |
---|
83 | // Copy constructor |
---|
84 | CbcHeuristicGreedyCover::CbcHeuristicGreedyCover(const CbcHeuristicGreedyCover & rhs) |
---|
85 | : |
---|
86 | CbcHeuristic(rhs), |
---|
87 | matrix_(rhs.matrix_), |
---|
88 | originalNumberRows_(rhs.originalNumberRows_), |
---|
89 | algorithm_(rhs.algorithm_), |
---|
90 | numberTimes_(rhs.numberTimes_) |
---|
91 | { |
---|
92 | } |
---|
93 | |
---|
94 | // Assignment operator |
---|
95 | CbcHeuristicGreedyCover & |
---|
96 | CbcHeuristicGreedyCover::operator=( const CbcHeuristicGreedyCover & rhs) |
---|
97 | { |
---|
98 | if (this != &rhs) { |
---|
99 | CbcHeuristic::operator=(rhs); |
---|
100 | matrix_ = rhs.matrix_; |
---|
101 | originalNumberRows_ = rhs.originalNumberRows_; |
---|
102 | algorithm_ = rhs.algorithm_; |
---|
103 | numberTimes_ = rhs.numberTimes_; |
---|
104 | } |
---|
105 | return *this; |
---|
106 | } |
---|
107 | // Returns 1 if solution, 0 if not |
---|
108 | int |
---|
109 | CbcHeuristicGreedyCover::solution(double & solutionValue, |
---|
110 | double * betterSolution) |
---|
111 | { |
---|
112 | numCouldRun_++; |
---|
113 | if (!model_) |
---|
114 | return 0; |
---|
115 | // See if to do |
---|
116 | if (!when() || (when() == 1 && model_->phase() != 1)) |
---|
117 | return 0; // switched off |
---|
118 | if (model_->getNodeCount() > numberTimes_) |
---|
119 | return 0; |
---|
120 | #ifdef HEURISTIC_INFORM |
---|
121 | printf("Entering heuristic %s - nRuns %d numCould %d when %d\n", |
---|
122 | heuristicName(),numRuns_,numCouldRun_,when_); |
---|
123 | #endif |
---|
124 | // See if at root node |
---|
125 | bool atRoot = model_->getNodeCount() == 0; |
---|
126 | int passNumber = model_->getCurrentPassNumber(); |
---|
127 | if (atRoot && passNumber > 1) |
---|
128 | return 0; |
---|
129 | OsiSolverInterface * solver = model_->solver(); |
---|
130 | const double * columnLower = solver->getColLower(); |
---|
131 | const double * columnUpper = solver->getColUpper(); |
---|
132 | // And original upper bounds in case we want to use them |
---|
133 | const double * originalUpper = model_->continuousSolver()->getColUpper(); |
---|
134 | // But not if algorithm says so |
---|
135 | if ((algorithm_ % 10) == 0) |
---|
136 | originalUpper = columnUpper; |
---|
137 | const double * rowLower = solver->getRowLower(); |
---|
138 | const double * solution = solver->getColSolution(); |
---|
139 | const double * objective = solver->getObjCoefficients(); |
---|
140 | double integerTolerance = model_->getDblParam(CbcModel::CbcIntegerTolerance); |
---|
141 | double primalTolerance; |
---|
142 | solver->getDblParam(OsiPrimalTolerance, primalTolerance); |
---|
143 | |
---|
144 | // This is number of rows when matrix was passed in |
---|
145 | int numberRows = originalNumberRows_; |
---|
146 | if (!numberRows) |
---|
147 | return 0; // switched off |
---|
148 | |
---|
149 | numRuns_++; |
---|
150 | assert (numberRows == matrix_.getNumRows()); |
---|
151 | int iRow, iColumn; |
---|
152 | double direction = solver->getObjSense(); |
---|
153 | double offset; |
---|
154 | solver->getDblParam(OsiObjOffset, offset); |
---|
155 | double newSolutionValue = -offset; |
---|
156 | int returnCode = 0; |
---|
157 | |
---|
158 | // Column copy |
---|
159 | const double * element = matrix_.getElements(); |
---|
160 | const int * row = matrix_.getIndices(); |
---|
161 | const CoinBigIndex * columnStart = matrix_.getVectorStarts(); |
---|
162 | const int * columnLength = matrix_.getVectorLengths(); |
---|
163 | |
---|
164 | // Get solution array for heuristic solution |
---|
165 | int numberColumns = solver->getNumCols(); |
---|
166 | double * newSolution = new double [numberColumns]; |
---|
167 | double * rowActivity = new double[numberRows]; |
---|
168 | memset(rowActivity, 0, numberRows*sizeof(double)); |
---|
169 | bool allOnes = true; |
---|
170 | // Get rounded down solution |
---|
171 | for (iColumn = 0; iColumn < numberColumns; iColumn++) { |
---|
172 | CoinBigIndex j; |
---|
173 | double value = solution[iColumn]; |
---|
174 | if (solver->isInteger(iColumn)) { |
---|
175 | // Round down integer |
---|
176 | if (fabs(floor(value + 0.5) - value) < integerTolerance) { |
---|
177 | value = floor(CoinMax(value + 1.0e-3, columnLower[iColumn])); |
---|
178 | } else { |
---|
179 | value = CoinMax(floor(value), columnLower[iColumn]); |
---|
180 | } |
---|
181 | } |
---|
182 | // make sure clean |
---|
183 | value = CoinMin(value, columnUpper[iColumn]); |
---|
184 | value = CoinMax(value, columnLower[iColumn]); |
---|
185 | newSolution[iColumn] = value; |
---|
186 | double cost = direction * objective[iColumn]; |
---|
187 | newSolutionValue += value * cost; |
---|
188 | for (j = columnStart[iColumn]; |
---|
189 | j < columnStart[iColumn] + columnLength[iColumn]; j++) { |
---|
190 | int iRow = row[j]; |
---|
191 | rowActivity[iRow] += value * element[j]; |
---|
192 | if (element[j] != 1.0) |
---|
193 | allOnes = false; |
---|
194 | } |
---|
195 | } |
---|
196 | // See if we round up |
---|
197 | bool roundup = ((algorithm_ % 100) != 0); |
---|
198 | if (roundup && allOnes) { |
---|
199 | // Get rounded up solution |
---|
200 | for (iColumn = 0; iColumn < numberColumns; iColumn++) { |
---|
201 | CoinBigIndex j; |
---|
202 | double value = solution[iColumn]; |
---|
203 | if (solver->isInteger(iColumn)) { |
---|
204 | // but round up if no activity |
---|
205 | if (roundup && value >= 0.499999 && !newSolution[iColumn]) { |
---|
206 | bool choose = true; |
---|
207 | for (j = columnStart[iColumn]; |
---|
208 | j < columnStart[iColumn] + columnLength[iColumn]; j++) { |
---|
209 | int iRow = row[j]; |
---|
210 | if (rowActivity[iRow]) { |
---|
211 | choose = false; |
---|
212 | break; |
---|
213 | } |
---|
214 | } |
---|
215 | if (choose) { |
---|
216 | newSolution[iColumn] = 1.0; |
---|
217 | double cost = direction * objective[iColumn]; |
---|
218 | newSolutionValue += cost; |
---|
219 | for (j = columnStart[iColumn]; |
---|
220 | j < columnStart[iColumn] + columnLength[iColumn]; j++) { |
---|
221 | int iRow = row[j]; |
---|
222 | rowActivity[iRow] += 1.0; |
---|
223 | } |
---|
224 | } |
---|
225 | } |
---|
226 | } |
---|
227 | } |
---|
228 | } |
---|
229 | // Get initial list |
---|
230 | int * which = new int [numberColumns]; |
---|
231 | for (iColumn = 0; iColumn < numberColumns; iColumn++) |
---|
232 | which[iColumn] = iColumn; |
---|
233 | int numberLook = numberColumns; |
---|
234 | // See if we want to perturb more |
---|
235 | double perturb = ((algorithm_ % 10) == 0) ? 0.1 : 0.25; |
---|
236 | // Keep going round until a solution |
---|
237 | while (true) { |
---|
238 | // Get column with best ratio |
---|
239 | int bestColumn = -1; |
---|
240 | double bestRatio = COIN_DBL_MAX; |
---|
241 | double bestStepSize = 0.0; |
---|
242 | int newNumber = 0; |
---|
243 | for (int jColumn = 0; jColumn < numberLook; jColumn++) { |
---|
244 | int iColumn = which[jColumn]; |
---|
245 | CoinBigIndex j; |
---|
246 | double value = newSolution[iColumn]; |
---|
247 | double cost = direction * objective[iColumn]; |
---|
248 | if (solver->isInteger(iColumn)) { |
---|
249 | // use current upper or original upper |
---|
250 | if (value + 0.99 < originalUpper[iColumn]) { |
---|
251 | double sum = 0.0; |
---|
252 | int numberExact = 0; |
---|
253 | for (j = columnStart[iColumn]; |
---|
254 | j < columnStart[iColumn] + columnLength[iColumn]; j++) { |
---|
255 | int iRow = row[j]; |
---|
256 | double gap = rowLower[iRow] - rowActivity[iRow]; |
---|
257 | double elementValue = allOnes ? 1.0 : element[j]; |
---|
258 | if (gap > 1.0e-7) { |
---|
259 | sum += CoinMin(elementValue, gap); |
---|
260 | if (fabs(elementValue - gap) < 1.0e-7) |
---|
261 | numberExact++; |
---|
262 | } |
---|
263 | } |
---|
264 | // could bias if exact |
---|
265 | if (sum > 0.0) { |
---|
266 | // add to next time |
---|
267 | which[newNumber++] = iColumn; |
---|
268 | double ratio = (cost / sum) * (1.0 + perturb * randomNumberGenerator_.randomDouble()); |
---|
269 | // If at root choose first |
---|
270 | if (atRoot) |
---|
271 | ratio = iColumn; |
---|
272 | if (ratio < bestRatio) { |
---|
273 | bestRatio = ratio; |
---|
274 | bestColumn = iColumn; |
---|
275 | bestStepSize = 1.0; |
---|
276 | } |
---|
277 | } |
---|
278 | } |
---|
279 | } else { |
---|
280 | // continuous |
---|
281 | if (value < columnUpper[iColumn]) { |
---|
282 | // Go through twice - first to get step length |
---|
283 | double step = 1.0e50; |
---|
284 | for (j = columnStart[iColumn]; |
---|
285 | j < columnStart[iColumn] + columnLength[iColumn]; j++) { |
---|
286 | int iRow = row[j]; |
---|
287 | if (rowActivity[iRow] < rowLower[iRow] - 1.0e-10 && |
---|
288 | element[j]*step + rowActivity[iRow] >= rowLower[iRow]) { |
---|
289 | step = (rowLower[iRow] - rowActivity[iRow]) / element[j];; |
---|
290 | } |
---|
291 | } |
---|
292 | // now ratio |
---|
293 | if (step < 1.0e50) { |
---|
294 | // add to next time |
---|
295 | which[newNumber++] = iColumn; |
---|
296 | assert (step > 0.0); |
---|
297 | double sum = 0.0; |
---|
298 | for (j = columnStart[iColumn]; |
---|
299 | j < columnStart[iColumn] + columnLength[iColumn]; j++) { |
---|
300 | int iRow = row[j]; |
---|
301 | double newActivity = element[j] * step + rowActivity[iRow]; |
---|
302 | if (rowActivity[iRow] < rowLower[iRow] - 1.0e-10 && |
---|
303 | newActivity >= rowLower[iRow] - 1.0e-12) { |
---|
304 | sum += element[j]; |
---|
305 | } |
---|
306 | } |
---|
307 | assert (sum > 0.0); |
---|
308 | double ratio = (cost / sum) * (1.0 + perturb * randomNumberGenerator_.randomDouble()); |
---|
309 | if (ratio < bestRatio) { |
---|
310 | bestRatio = ratio; |
---|
311 | bestColumn = iColumn; |
---|
312 | bestStepSize = step; |
---|
313 | } |
---|
314 | } |
---|
315 | } |
---|
316 | } |
---|
317 | } |
---|
318 | if (bestColumn < 0) |
---|
319 | break; // we have finished |
---|
320 | // Increase chosen column |
---|
321 | newSolution[bestColumn] += bestStepSize; |
---|
322 | double cost = direction * objective[bestColumn]; |
---|
323 | newSolutionValue += bestStepSize * cost; |
---|
324 | for (CoinBigIndex j = columnStart[bestColumn]; |
---|
325 | j < columnStart[bestColumn] + columnLength[bestColumn]; j++) { |
---|
326 | int iRow = row[j]; |
---|
327 | rowActivity[iRow] += bestStepSize * element[j]; |
---|
328 | } |
---|
329 | } |
---|
330 | delete [] which; |
---|
331 | if (newSolutionValue < solutionValue) { |
---|
332 | // check feasible |
---|
333 | memset(rowActivity, 0, numberRows*sizeof(double)); |
---|
334 | for (iColumn = 0; iColumn < numberColumns; iColumn++) { |
---|
335 | CoinBigIndex j; |
---|
336 | double value = newSolution[iColumn]; |
---|
337 | if (value) { |
---|
338 | for (j = columnStart[iColumn]; |
---|
339 | j < columnStart[iColumn] + columnLength[iColumn]; j++) { |
---|
340 | int iRow = row[j]; |
---|
341 | rowActivity[iRow] += value * element[j]; |
---|
342 | } |
---|
343 | } |
---|
344 | } |
---|
345 | // check was approximately feasible |
---|
346 | bool feasible = true; |
---|
347 | for (iRow = 0; iRow < numberRows; iRow++) { |
---|
348 | if (rowActivity[iRow] < rowLower[iRow]) { |
---|
349 | if (rowActivity[iRow] < rowLower[iRow] - 10.0*primalTolerance) |
---|
350 | feasible = false; |
---|
351 | } |
---|
352 | } |
---|
353 | if (feasible) { |
---|
354 | // new solution |
---|
355 | memcpy(betterSolution, newSolution, numberColumns*sizeof(double)); |
---|
356 | solutionValue = newSolutionValue; |
---|
357 | //printf("** Solution of %g found by rounding\n",newSolutionValue); |
---|
358 | returnCode = 1; |
---|
359 | } else { |
---|
360 | // Can easily happen |
---|
361 | //printf("Debug CbcHeuristicGreedyCover giving bad solution\n"); |
---|
362 | } |
---|
363 | } |
---|
364 | delete [] newSolution; |
---|
365 | delete [] rowActivity; |
---|
366 | return returnCode; |
---|
367 | } |
---|
368 | // update model |
---|
369 | void CbcHeuristicGreedyCover::setModel(CbcModel * model) |
---|
370 | { |
---|
371 | gutsOfConstructor(model); |
---|
372 | validate(); |
---|
373 | } |
---|
374 | // Resets stuff if model changes |
---|
375 | void |
---|
376 | CbcHeuristicGreedyCover::resetModel(CbcModel * model) |
---|
377 | { |
---|
378 | gutsOfConstructor(model); |
---|
379 | } |
---|
380 | // Validate model i.e. sets when_ to 0 if necessary (may be NULL) |
---|
381 | void |
---|
382 | CbcHeuristicGreedyCover::validate() |
---|
383 | { |
---|
384 | if (model_ && when() < 10) { |
---|
385 | if (model_->numberIntegers() != |
---|
386 | model_->numberObjects() && (model_->numberObjects() || |
---|
387 | (model_->specialOptions()&1024) == 0)) { |
---|
388 | int numberOdd = 0; |
---|
389 | for (int i = 0; i < model_->numberObjects(); i++) { |
---|
390 | if (!model_->object(i)->canDoHeuristics()) |
---|
391 | numberOdd++; |
---|
392 | } |
---|
393 | if (numberOdd) |
---|
394 | setWhen(0); |
---|
395 | } |
---|
396 | // Only works if costs positive, coefficients positive and all rows G |
---|
397 | OsiSolverInterface * solver = model_->solver(); |
---|
398 | const double * columnLower = solver->getColLower(); |
---|
399 | const double * rowUpper = solver->getRowUpper(); |
---|
400 | const double * objective = solver->getObjCoefficients(); |
---|
401 | double direction = solver->getObjSense(); |
---|
402 | |
---|
403 | int numberRows = solver->getNumRows(); |
---|
404 | int numberColumns = solver->getNumCols(); |
---|
405 | // Column copy |
---|
406 | matrix_.setDimensions(numberRows,numberColumns); |
---|
407 | const double * element = matrix_.getElements(); |
---|
408 | const CoinBigIndex * columnStart = matrix_.getVectorStarts(); |
---|
409 | const int * columnLength = matrix_.getVectorLengths(); |
---|
410 | bool good = true; |
---|
411 | for (int iRow = 0; iRow < numberRows; iRow++) { |
---|
412 | if (rowUpper[iRow] < 1.0e30) |
---|
413 | good = false; |
---|
414 | } |
---|
415 | for (int iColumn = 0; iColumn < numberColumns; iColumn++) { |
---|
416 | if (objective[iColumn]*direction < 0.0) |
---|
417 | good = false; |
---|
418 | if (columnLower[iColumn] < 0.0) |
---|
419 | good = false; |
---|
420 | CoinBigIndex j; |
---|
421 | for (j = columnStart[iColumn]; |
---|
422 | j < columnStart[iColumn] + columnLength[iColumn]; j++) { |
---|
423 | if (element[j] < 0.0) |
---|
424 | good = false; |
---|
425 | } |
---|
426 | } |
---|
427 | if (!good) |
---|
428 | setWhen(0); // switch off |
---|
429 | } |
---|
430 | } |
---|
431 | // Default Constructor |
---|
432 | CbcHeuristicGreedyEquality::CbcHeuristicGreedyEquality() |
---|
433 | : CbcHeuristic() |
---|
434 | { |
---|
435 | // matrix will automatically be empty |
---|
436 | fraction_ = 1.0; // no branch and bound |
---|
437 | originalNumberRows_ = 0; |
---|
438 | algorithm_ = 0; |
---|
439 | numberTimes_ = 100; |
---|
440 | whereFrom_ = 1; |
---|
441 | } |
---|
442 | |
---|
443 | // Constructor from model |
---|
444 | CbcHeuristicGreedyEquality::CbcHeuristicGreedyEquality(CbcModel & model) |
---|
445 | : CbcHeuristic(model) |
---|
446 | { |
---|
447 | // Get a copy of original matrix |
---|
448 | gutsOfConstructor(&model); |
---|
449 | fraction_ = 1.0; // no branch and bound |
---|
450 | algorithm_ = 0; |
---|
451 | numberTimes_ = 100; |
---|
452 | whereFrom_ = 1; |
---|
453 | } |
---|
454 | |
---|
455 | // Destructor |
---|
456 | CbcHeuristicGreedyEquality::~CbcHeuristicGreedyEquality () |
---|
457 | { |
---|
458 | } |
---|
459 | |
---|
460 | // Clone |
---|
461 | CbcHeuristic * |
---|
462 | CbcHeuristicGreedyEquality::clone() const |
---|
463 | { |
---|
464 | return new CbcHeuristicGreedyEquality(*this); |
---|
465 | } |
---|
466 | // Guts of constructor from a CbcModel |
---|
467 | void |
---|
468 | CbcHeuristicGreedyEquality::gutsOfConstructor(CbcModel * model) |
---|
469 | { |
---|
470 | model_ = model; |
---|
471 | // Get a copy of original matrix |
---|
472 | assert(model->solver()); |
---|
473 | if (model->solver()->getNumRows()) { |
---|
474 | matrix_ = *model->solver()->getMatrixByCol(); |
---|
475 | } |
---|
476 | originalNumberRows_ = model->solver()->getNumRows(); |
---|
477 | } |
---|
478 | // Create C++ lines to get to current state |
---|
479 | void |
---|
480 | CbcHeuristicGreedyEquality::generateCpp( FILE * fp) |
---|
481 | { |
---|
482 | CbcHeuristicGreedyEquality other; |
---|
483 | fprintf(fp, "0#include \"CbcHeuristicGreedy.hpp\"\n"); |
---|
484 | fprintf(fp, "3 CbcHeuristicGreedyEquality heuristicGreedyEquality(*cbcModel);\n"); |
---|
485 | CbcHeuristic::generateCpp(fp, "heuristicGreedyEquality"); |
---|
486 | if (algorithm_ != other.algorithm_) |
---|
487 | fprintf(fp, "3 heuristicGreedyEquality.setAlgorithm(%d);\n", algorithm_); |
---|
488 | else |
---|
489 | fprintf(fp, "4 heuristicGreedyEquality.setAlgorithm(%d);\n", algorithm_); |
---|
490 | if (fraction_ != other.fraction_) |
---|
491 | fprintf(fp, "3 heuristicGreedyEquality.setFraction(%g);\n", fraction_); |
---|
492 | else |
---|
493 | fprintf(fp, "4 heuristicGreedyEquality.setFraction(%g);\n", fraction_); |
---|
494 | if (numberTimes_ != other.numberTimes_) |
---|
495 | fprintf(fp, "3 heuristicGreedyEquality.setNumberTimes(%d);\n", numberTimes_); |
---|
496 | else |
---|
497 | fprintf(fp, "4 heuristicGreedyEquality.setNumberTimes(%d);\n", numberTimes_); |
---|
498 | fprintf(fp, "3 cbcModel->addHeuristic(&heuristicGreedyEquality);\n"); |
---|
499 | } |
---|
500 | |
---|
501 | // Copy constructor |
---|
502 | CbcHeuristicGreedyEquality::CbcHeuristicGreedyEquality(const CbcHeuristicGreedyEquality & rhs) |
---|
503 | : |
---|
504 | CbcHeuristic(rhs), |
---|
505 | matrix_(rhs.matrix_), |
---|
506 | fraction_(rhs.fraction_), |
---|
507 | originalNumberRows_(rhs.originalNumberRows_), |
---|
508 | algorithm_(rhs.algorithm_), |
---|
509 | numberTimes_(rhs.numberTimes_) |
---|
510 | { |
---|
511 | } |
---|
512 | |
---|
513 | // Assignment operator |
---|
514 | CbcHeuristicGreedyEquality & |
---|
515 | CbcHeuristicGreedyEquality::operator=( const CbcHeuristicGreedyEquality & rhs) |
---|
516 | { |
---|
517 | if (this != &rhs) { |
---|
518 | CbcHeuristic::operator=(rhs); |
---|
519 | matrix_ = rhs.matrix_; |
---|
520 | fraction_ = rhs.fraction_; |
---|
521 | originalNumberRows_ = rhs.originalNumberRows_; |
---|
522 | algorithm_ = rhs.algorithm_; |
---|
523 | numberTimes_ = rhs.numberTimes_; |
---|
524 | } |
---|
525 | return *this; |
---|
526 | } |
---|
527 | // Returns 1 if solution, 0 if not |
---|
528 | int |
---|
529 | CbcHeuristicGreedyEquality::solution(double & solutionValue, |
---|
530 | double * betterSolution) |
---|
531 | { |
---|
532 | numCouldRun_++; |
---|
533 | if (!model_) |
---|
534 | return 0; |
---|
535 | // See if to do |
---|
536 | if (!when() || (when() == 1 && model_->phase() != 1)) |
---|
537 | return 0; // switched off |
---|
538 | if (model_->getNodeCount() > numberTimes_) |
---|
539 | return 0; |
---|
540 | // See if at root node |
---|
541 | bool atRoot = model_->getNodeCount() == 0; |
---|
542 | int passNumber = model_->getCurrentPassNumber(); |
---|
543 | if (atRoot && passNumber > 1) |
---|
544 | return 0; |
---|
545 | OsiSolverInterface * solver = model_->solver(); |
---|
546 | const double * columnLower = solver->getColLower(); |
---|
547 | const double * columnUpper = solver->getColUpper(); |
---|
548 | // And original upper bounds in case we want to use them |
---|
549 | const double * originalUpper = model_->continuousSolver()->getColUpper(); |
---|
550 | // But not if algorithm says so |
---|
551 | if ((algorithm_ % 10) == 0) |
---|
552 | originalUpper = columnUpper; |
---|
553 | const double * rowLower = solver->getRowLower(); |
---|
554 | const double * rowUpper = solver->getRowUpper(); |
---|
555 | const double * solution = solver->getColSolution(); |
---|
556 | const double * objective = solver->getObjCoefficients(); |
---|
557 | double integerTolerance = model_->getDblParam(CbcModel::CbcIntegerTolerance); |
---|
558 | double primalTolerance; |
---|
559 | solver->getDblParam(OsiPrimalTolerance, primalTolerance); |
---|
560 | |
---|
561 | // This is number of rows when matrix was passed in |
---|
562 | int numberRows = originalNumberRows_; |
---|
563 | if (!numberRows) |
---|
564 | return 0; // switched off |
---|
565 | numRuns_++; |
---|
566 | |
---|
567 | assert (numberRows == matrix_.getNumRows()); |
---|
568 | int iRow, iColumn; |
---|
569 | double direction = solver->getObjSense(); |
---|
570 | double offset; |
---|
571 | solver->getDblParam(OsiObjOffset, offset); |
---|
572 | double newSolutionValue = -offset; |
---|
573 | int returnCode = 0; |
---|
574 | |
---|
575 | // Column copy |
---|
576 | const double * element = matrix_.getElements(); |
---|
577 | const int * row = matrix_.getIndices(); |
---|
578 | const CoinBigIndex * columnStart = matrix_.getVectorStarts(); |
---|
579 | const int * columnLength = matrix_.getVectorLengths(); |
---|
580 | |
---|
581 | // Get solution array for heuristic solution |
---|
582 | int numberColumns = solver->getNumCols(); |
---|
583 | double * newSolution = new double [numberColumns]; |
---|
584 | double * rowActivity = new double[numberRows]; |
---|
585 | memset(rowActivity, 0, numberRows*sizeof(double)); |
---|
586 | double rhsNeeded = 0; |
---|
587 | for (iRow = 0; iRow < numberRows; iRow++) |
---|
588 | rhsNeeded += rowUpper[iRow]; |
---|
589 | rhsNeeded *= fraction_; |
---|
590 | bool allOnes = true; |
---|
591 | // Get rounded down solution |
---|
592 | for (iColumn = 0; iColumn < numberColumns; iColumn++) { |
---|
593 | CoinBigIndex j; |
---|
594 | double value = solution[iColumn]; |
---|
595 | if (solver->isInteger(iColumn)) { |
---|
596 | // Round down integer |
---|
597 | if (fabs(floor(value + 0.5) - value) < integerTolerance) { |
---|
598 | value = floor(CoinMax(value + 1.0e-3, columnLower[iColumn])); |
---|
599 | } else { |
---|
600 | value = CoinMax(floor(value), columnLower[iColumn]); |
---|
601 | } |
---|
602 | } |
---|
603 | // make sure clean |
---|
604 | value = CoinMin(value, columnUpper[iColumn]); |
---|
605 | value = CoinMax(value, columnLower[iColumn]); |
---|
606 | newSolution[iColumn] = value; |
---|
607 | double cost = direction * objective[iColumn]; |
---|
608 | newSolutionValue += value * cost; |
---|
609 | for (j = columnStart[iColumn]; |
---|
610 | j < columnStart[iColumn] + columnLength[iColumn]; j++) { |
---|
611 | int iRow = row[j]; |
---|
612 | rowActivity[iRow] += value * element[j]; |
---|
613 | rhsNeeded -= value * element[j]; |
---|
614 | if (element[j] != 1.0) |
---|
615 | allOnes = false; |
---|
616 | } |
---|
617 | } |
---|
618 | // See if we round up |
---|
619 | bool roundup = ((algorithm_ % 100) != 0); |
---|
620 | if (roundup && allOnes) { |
---|
621 | // Get rounded up solution |
---|
622 | for (iColumn = 0; iColumn < numberColumns; iColumn++) { |
---|
623 | CoinBigIndex j; |
---|
624 | double value = solution[iColumn]; |
---|
625 | if (solver->isInteger(iColumn)) { |
---|
626 | // but round up if no activity |
---|
627 | if (roundup && value >= 0.6 && !newSolution[iColumn]) { |
---|
628 | bool choose = true; |
---|
629 | for (j = columnStart[iColumn]; |
---|
630 | j < columnStart[iColumn] + columnLength[iColumn]; j++) { |
---|
631 | int iRow = row[j]; |
---|
632 | if (rowActivity[iRow]) { |
---|
633 | choose = false; |
---|
634 | break; |
---|
635 | } |
---|
636 | } |
---|
637 | if (choose) { |
---|
638 | newSolution[iColumn] = 1.0; |
---|
639 | double cost = direction * objective[iColumn]; |
---|
640 | newSolutionValue += cost; |
---|
641 | for (j = columnStart[iColumn]; |
---|
642 | j < columnStart[iColumn] + columnLength[iColumn]; j++) { |
---|
643 | int iRow = row[j]; |
---|
644 | rowActivity[iRow] += 1.0; |
---|
645 | rhsNeeded -= 1.0; |
---|
646 | } |
---|
647 | } |
---|
648 | } |
---|
649 | } |
---|
650 | } |
---|
651 | } |
---|
652 | // Get initial list |
---|
653 | int * which = new int [numberColumns]; |
---|
654 | for (iColumn = 0; iColumn < numberColumns; iColumn++) |
---|
655 | which[iColumn] = iColumn; |
---|
656 | int numberLook = numberColumns; |
---|
657 | // See if we want to perturb more |
---|
658 | double perturb = ((algorithm_ % 10) == 0) ? 0.1 : 0.25; |
---|
659 | // Keep going round until a solution |
---|
660 | while (true) { |
---|
661 | // Get column with best ratio |
---|
662 | int bestColumn = -1; |
---|
663 | double bestRatio = COIN_DBL_MAX; |
---|
664 | double bestStepSize = 0.0; |
---|
665 | int newNumber = 0; |
---|
666 | for (int jColumn = 0; jColumn < numberLook; jColumn++) { |
---|
667 | int iColumn = which[jColumn]; |
---|
668 | CoinBigIndex j; |
---|
669 | double value = newSolution[iColumn]; |
---|
670 | double cost = direction * objective[iColumn]; |
---|
671 | if (solver->isInteger(iColumn)) { |
---|
672 | // use current upper or original upper |
---|
673 | if (value + 0.9999 < originalUpper[iColumn]) { |
---|
674 | double movement = 1.0; |
---|
675 | double sum = 0.0; |
---|
676 | for (j = columnStart[iColumn]; |
---|
677 | j < columnStart[iColumn] + columnLength[iColumn]; j++) { |
---|
678 | int iRow = row[j]; |
---|
679 | double gap = rowUpper[iRow] - rowActivity[iRow]; |
---|
680 | double elementValue = allOnes ? 1.0 : element[j]; |
---|
681 | sum += elementValue; |
---|
682 | if (movement*elementValue > gap) { |
---|
683 | movement = gap / elementValue; |
---|
684 | } |
---|
685 | } |
---|
686 | if (movement > 0.999999) { |
---|
687 | // add to next time |
---|
688 | which[newNumber++] = iColumn; |
---|
689 | double ratio = (cost / sum) * (1.0 + perturb * randomNumberGenerator_.randomDouble()); |
---|
690 | // If at root |
---|
691 | if (atRoot) { |
---|
692 | if (fraction_ == 1.0) |
---|
693 | ratio = iColumn; // choose first |
---|
694 | else |
---|
695 | ratio = - solution[iColumn]; // choose largest |
---|
696 | } |
---|
697 | if (ratio < bestRatio) { |
---|
698 | bestRatio = ratio; |
---|
699 | bestColumn = iColumn; |
---|
700 | bestStepSize = 1.0; |
---|
701 | } |
---|
702 | } |
---|
703 | } |
---|
704 | } else { |
---|
705 | // continuous |
---|
706 | if (value < columnUpper[iColumn]) { |
---|
707 | double movement = 1.0e50; |
---|
708 | double sum = 0.0; |
---|
709 | for (j = columnStart[iColumn]; |
---|
710 | j < columnStart[iColumn] + columnLength[iColumn]; j++) { |
---|
711 | int iRow = row[j]; |
---|
712 | if (element[j]*movement + rowActivity[iRow] > rowUpper[iRow]) { |
---|
713 | movement = (rowUpper[iRow] - rowActivity[iRow]) / element[j];; |
---|
714 | } |
---|
715 | sum += element[j]; |
---|
716 | } |
---|
717 | // now ratio |
---|
718 | if (movement > 1.0e-7) { |
---|
719 | // add to next time |
---|
720 | which[newNumber++] = iColumn; |
---|
721 | double ratio = (cost / sum) * (1.0 + perturb * randomNumberGenerator_.randomDouble()); |
---|
722 | if (ratio < bestRatio) { |
---|
723 | bestRatio = ratio; |
---|
724 | bestColumn = iColumn; |
---|
725 | bestStepSize = movement; |
---|
726 | } |
---|
727 | } |
---|
728 | } |
---|
729 | } |
---|
730 | } |
---|
731 | if (bestColumn < 0) |
---|
732 | break; // we have finished |
---|
733 | // Increase chosen column |
---|
734 | newSolution[bestColumn] += bestStepSize; |
---|
735 | double cost = direction * objective[bestColumn]; |
---|
736 | newSolutionValue += bestStepSize * cost; |
---|
737 | for (CoinBigIndex j = columnStart[bestColumn]; |
---|
738 | j < columnStart[bestColumn] + columnLength[bestColumn]; j++) { |
---|
739 | int iRow = row[j]; |
---|
740 | rowActivity[iRow] += bestStepSize * element[j]; |
---|
741 | rhsNeeded -= bestStepSize * element[j]; |
---|
742 | } |
---|
743 | if (rhsNeeded < 1.0e-8) |
---|
744 | break; |
---|
745 | } |
---|
746 | delete [] which; |
---|
747 | if (fraction_ < 1.0 && rhsNeeded < 1.0e-8 && newSolutionValue < solutionValue) { |
---|
748 | // do branch and cut |
---|
749 | // fix all nonzero |
---|
750 | OsiSolverInterface * newSolver = model_->continuousSolver()->clone(); |
---|
751 | for (iColumn = 0; iColumn < numberColumns; iColumn++) { |
---|
752 | if (newSolver->isInteger(iColumn)) |
---|
753 | newSolver->setColLower(iColumn, newSolution[iColumn]); |
---|
754 | } |
---|
755 | int returnCode = smallBranchAndBound(newSolver, 200, newSolution, newSolutionValue, |
---|
756 | solutionValue, "CbcHeuristicGreedy"); |
---|
757 | if (returnCode < 0) |
---|
758 | returnCode = 0; // returned on size |
---|
759 | if ((returnCode&2) != 0) { |
---|
760 | // could add cut |
---|
761 | returnCode &= ~2; |
---|
762 | } |
---|
763 | rhsNeeded = 1.0 - returnCode; |
---|
764 | delete newSolver; |
---|
765 | } |
---|
766 | if (newSolutionValue < solutionValue && rhsNeeded < 1.0e-8) { |
---|
767 | // check feasible |
---|
768 | memset(rowActivity, 0, numberRows*sizeof(double)); |
---|
769 | for (iColumn = 0; iColumn < numberColumns; iColumn++) { |
---|
770 | CoinBigIndex j; |
---|
771 | double value = newSolution[iColumn]; |
---|
772 | if (value) { |
---|
773 | for (j = columnStart[iColumn]; |
---|
774 | j < columnStart[iColumn] + columnLength[iColumn]; j++) { |
---|
775 | int iRow = row[j]; |
---|
776 | rowActivity[iRow] += value * element[j]; |
---|
777 | } |
---|
778 | } |
---|
779 | } |
---|
780 | // check was approximately feasible |
---|
781 | bool feasible = true; |
---|
782 | for (iRow = 0; iRow < numberRows; iRow++) { |
---|
783 | if (rowActivity[iRow] < rowLower[iRow]) { |
---|
784 | if (rowActivity[iRow] < rowLower[iRow] - 10.0*primalTolerance) |
---|
785 | feasible = false; |
---|
786 | } |
---|
787 | } |
---|
788 | if (feasible) { |
---|
789 | // new solution |
---|
790 | memcpy(betterSolution, newSolution, numberColumns*sizeof(double)); |
---|
791 | solutionValue = newSolutionValue; |
---|
792 | returnCode = 1; |
---|
793 | } |
---|
794 | } |
---|
795 | delete [] newSolution; |
---|
796 | delete [] rowActivity; |
---|
797 | if (atRoot && fraction_ == 1.0) { |
---|
798 | // try quick search |
---|
799 | fraction_ = 0.4; |
---|
800 | int newCode = this->solution(solutionValue, betterSolution); |
---|
801 | if (newCode) |
---|
802 | returnCode = 1; |
---|
803 | fraction_ = 1.0; |
---|
804 | } |
---|
805 | return returnCode; |
---|
806 | } |
---|
807 | // update model |
---|
808 | void CbcHeuristicGreedyEquality::setModel(CbcModel * model) |
---|
809 | { |
---|
810 | gutsOfConstructor(model); |
---|
811 | validate(); |
---|
812 | } |
---|
813 | // Resets stuff if model changes |
---|
814 | void |
---|
815 | CbcHeuristicGreedyEquality::resetModel(CbcModel * model) |
---|
816 | { |
---|
817 | gutsOfConstructor(model); |
---|
818 | } |
---|
819 | // Validate model i.e. sets when_ to 0 if necessary (may be NULL) |
---|
820 | void |
---|
821 | CbcHeuristicGreedyEquality::validate() |
---|
822 | { |
---|
823 | if (model_ && when() < 10) { |
---|
824 | if (model_->numberIntegers() != |
---|
825 | model_->numberObjects()) |
---|
826 | setWhen(0); |
---|
827 | // Only works if costs positive, coefficients positive and all rows E or L |
---|
828 | // And if values are integer |
---|
829 | OsiSolverInterface * solver = model_->solver(); |
---|
830 | const double * columnLower = solver->getColLower(); |
---|
831 | const double * rowUpper = solver->getRowUpper(); |
---|
832 | const double * rowLower = solver->getRowLower(); |
---|
833 | const double * objective = solver->getObjCoefficients(); |
---|
834 | double direction = solver->getObjSense(); |
---|
835 | |
---|
836 | int numberRows = solver->getNumRows(); |
---|
837 | int numberColumns = solver->getNumCols(); |
---|
838 | matrix_.setDimensions(numberRows,numberColumns); |
---|
839 | // Column copy |
---|
840 | const double * element = matrix_.getElements(); |
---|
841 | const CoinBigIndex * columnStart = matrix_.getVectorStarts(); |
---|
842 | const int * columnLength = matrix_.getVectorLengths(); |
---|
843 | bool good = true; |
---|
844 | for (int iRow = 0; iRow < numberRows; iRow++) { |
---|
845 | if (rowUpper[iRow] > 1.0e30) |
---|
846 | good = false; |
---|
847 | if (rowLower[iRow] > 0.0 && rowLower[iRow] != rowUpper[iRow]) |
---|
848 | good = false; |
---|
849 | if (floor(rowUpper[iRow] + 0.5) != rowUpper[iRow]) |
---|
850 | good = false; |
---|
851 | } |
---|
852 | for (int iColumn = 0; iColumn < numberColumns; iColumn++) { |
---|
853 | if (objective[iColumn]*direction < 0.0) |
---|
854 | good = false; |
---|
855 | if (columnLower[iColumn] < 0.0) |
---|
856 | good = false; |
---|
857 | CoinBigIndex j; |
---|
858 | for (j = columnStart[iColumn]; |
---|
859 | j < columnStart[iColumn] + columnLength[iColumn]; j++) { |
---|
860 | if (element[j] < 0.0) |
---|
861 | good = false; |
---|
862 | if (floor(element[j] + 0.5) != element[j]) |
---|
863 | good = false; |
---|
864 | } |
---|
865 | } |
---|
866 | if (!good) |
---|
867 | setWhen(0); // switch off |
---|
868 | } |
---|
869 | } |
---|
870 | |
---|
871 | // Default Constructor |
---|
872 | CbcHeuristicGreedySOS::CbcHeuristicGreedySOS() |
---|
873 | : CbcHeuristic() |
---|
874 | { |
---|
875 | originalRhs_ = NULL; |
---|
876 | // matrix will automatically be empty |
---|
877 | originalNumberRows_ = 0; |
---|
878 | algorithm_ = 0; |
---|
879 | numberTimes_ = 100; |
---|
880 | } |
---|
881 | |
---|
882 | // Constructor from model |
---|
883 | CbcHeuristicGreedySOS::CbcHeuristicGreedySOS(CbcModel & model) |
---|
884 | : CbcHeuristic(model) |
---|
885 | { |
---|
886 | gutsOfConstructor(&model); |
---|
887 | algorithm_ = 2; |
---|
888 | numberTimes_ = 100; |
---|
889 | whereFrom_ = 1; |
---|
890 | } |
---|
891 | |
---|
892 | // Destructor |
---|
893 | CbcHeuristicGreedySOS::~CbcHeuristicGreedySOS () |
---|
894 | { |
---|
895 | delete [] originalRhs_; |
---|
896 | } |
---|
897 | |
---|
898 | // Clone |
---|
899 | CbcHeuristic * |
---|
900 | CbcHeuristicGreedySOS::clone() const |
---|
901 | { |
---|
902 | return new CbcHeuristicGreedySOS(*this); |
---|
903 | } |
---|
904 | // Guts of constructor from a CbcModel |
---|
905 | void |
---|
906 | CbcHeuristicGreedySOS::gutsOfConstructor(CbcModel * model) |
---|
907 | { |
---|
908 | model_ = model; |
---|
909 | // Get a copy of original matrix |
---|
910 | assert(model->solver()); |
---|
911 | if (model->solver()->getNumRows()) { |
---|
912 | matrix_ = *model->solver()->getMatrixByCol(); |
---|
913 | } |
---|
914 | originalNumberRows_ = model->solver()->getNumRows(); |
---|
915 | originalRhs_ = new double [originalNumberRows_]; |
---|
916 | } |
---|
917 | // Create C++ lines to get to current state |
---|
918 | void |
---|
919 | CbcHeuristicGreedySOS::generateCpp( FILE * fp) |
---|
920 | { |
---|
921 | CbcHeuristicGreedySOS other; |
---|
922 | fprintf(fp, "0#include \"CbcHeuristicGreedy.hpp\"\n"); |
---|
923 | fprintf(fp, "3 CbcHeuristicGreedySOS heuristicGreedySOS(*cbcModel);\n"); |
---|
924 | CbcHeuristic::generateCpp(fp, "heuristicGreedySOS"); |
---|
925 | if (algorithm_ != other.algorithm_) |
---|
926 | fprintf(fp, "3 heuristicGreedySOS.setAlgorithm(%d);\n", algorithm_); |
---|
927 | else |
---|
928 | fprintf(fp, "4 heuristicGreedySOS.setAlgorithm(%d);\n", algorithm_); |
---|
929 | if (numberTimes_ != other.numberTimes_) |
---|
930 | fprintf(fp, "3 heuristicGreedySOS.setNumberTimes(%d);\n", numberTimes_); |
---|
931 | else |
---|
932 | fprintf(fp, "4 heuristicGreedySOS.setNumberTimes(%d);\n", numberTimes_); |
---|
933 | fprintf(fp, "3 cbcModel->addHeuristic(&heuristicGreedySOS);\n"); |
---|
934 | } |
---|
935 | |
---|
936 | // Copy constructor |
---|
937 | CbcHeuristicGreedySOS::CbcHeuristicGreedySOS(const CbcHeuristicGreedySOS & rhs) |
---|
938 | : |
---|
939 | CbcHeuristic(rhs), |
---|
940 | matrix_(rhs.matrix_), |
---|
941 | originalNumberRows_(rhs.originalNumberRows_), |
---|
942 | algorithm_(rhs.algorithm_), |
---|
943 | numberTimes_(rhs.numberTimes_) |
---|
944 | { |
---|
945 | originalRhs_ = CoinCopyOfArray(rhs.originalRhs_,originalNumberRows_); |
---|
946 | } |
---|
947 | |
---|
948 | // Assignment operator |
---|
949 | CbcHeuristicGreedySOS & |
---|
950 | CbcHeuristicGreedySOS::operator=( const CbcHeuristicGreedySOS & rhs) |
---|
951 | { |
---|
952 | if (this != &rhs) { |
---|
953 | CbcHeuristic::operator=(rhs); |
---|
954 | matrix_ = rhs.matrix_; |
---|
955 | originalNumberRows_ = rhs.originalNumberRows_; |
---|
956 | algorithm_ = rhs.algorithm_; |
---|
957 | numberTimes_ = rhs.numberTimes_; |
---|
958 | delete [] originalRhs_; |
---|
959 | originalRhs_ = CoinCopyOfArray(rhs.originalRhs_,originalNumberRows_); |
---|
960 | } |
---|
961 | return *this; |
---|
962 | } |
---|
963 | // Returns 1 if solution, 0 if not |
---|
964 | int |
---|
965 | CbcHeuristicGreedySOS::solution(double & solutionValue, |
---|
966 | double * betterSolution) |
---|
967 | { |
---|
968 | numCouldRun_++; |
---|
969 | if (!model_) |
---|
970 | return 0; |
---|
971 | // See if to do |
---|
972 | if (!when() || (when() == 1 && model_->phase() != 1)) |
---|
973 | return 0; // switched off |
---|
974 | if (model_->getNodeCount() > numberTimes_) |
---|
975 | return 0; |
---|
976 | // See if at root node |
---|
977 | bool atRoot = model_->getNodeCount() == 0; |
---|
978 | int passNumber = model_->getCurrentPassNumber(); |
---|
979 | if (atRoot && passNumber > 1) |
---|
980 | return 0; |
---|
981 | OsiSolverInterface * solver = model_->solver(); |
---|
982 | int numberColumns = solver->getNumCols(); |
---|
983 | // This is number of rows when matrix was passed in |
---|
984 | int numberRows = originalNumberRows_; |
---|
985 | if (!numberRows) |
---|
986 | return 0; // switched off |
---|
987 | |
---|
988 | const double * columnLower = solver->getColLower(); |
---|
989 | const double * columnUpper = solver->getColUpper(); |
---|
990 | // modified rhs |
---|
991 | double * rhs = CoinCopyOfArray(originalRhs_,numberRows); |
---|
992 | // Column copy |
---|
993 | const double * element = matrix_.getElements(); |
---|
994 | const int * row = matrix_.getIndices(); |
---|
995 | const CoinBigIndex * columnStart = matrix_.getVectorStarts(); |
---|
996 | const int * columnLength = matrix_.getVectorLengths(); |
---|
997 | int * sosRow = new int [numberColumns]; |
---|
998 | char * sos = new char [numberRows]; |
---|
999 | memset(sos,'a',numberRows); |
---|
1000 | int nonSOS=0; |
---|
1001 | // If bit set then use current |
---|
1002 | if ((algorithm_&1)!=0) { |
---|
1003 | const CoinPackedMatrix * matrix = solver->getMatrixByCol(); |
---|
1004 | element = matrix->getElements(); |
---|
1005 | row = matrix->getIndices(); |
---|
1006 | columnStart = matrix->getVectorStarts(); |
---|
1007 | columnLength = matrix->getVectorLengths(); |
---|
1008 | //rhs = new double [numberRows]; |
---|
1009 | const double * rowLower = solver->getRowLower(); |
---|
1010 | const double * rowUpper = solver->getRowUpper(); |
---|
1011 | bool good = true; |
---|
1012 | for (int iRow = 0; iRow < numberRows; iRow++) { |
---|
1013 | if (rowLower[iRow] == 1.0 && rowUpper[iRow] == 1.0) { |
---|
1014 | // SOS |
---|
1015 | rhs[iRow]=-1.0; |
---|
1016 | } else if (rowLower[iRow] > 0.0 && rowUpper[iRow] < 1.0e10) { |
---|
1017 | good = false; |
---|
1018 | } else if (rowUpper[iRow] < 0.0) { |
---|
1019 | good = false; |
---|
1020 | } else if (rowUpper[iRow] < 1.0e10) { |
---|
1021 | rhs[iRow]=rowUpper[iRow]; |
---|
1022 | if (rhs[iRow]<0) |
---|
1023 | sos[iRow]=0; // can't be SOS |
---|
1024 | } else { |
---|
1025 | rhs[iRow]=rowLower[iRow]; |
---|
1026 | if (rhs[iRow]<0) |
---|
1027 | sos[iRow]=0; // can't be SOS |
---|
1028 | } |
---|
1029 | } |
---|
1030 | for (int iColumn = 0; iColumn < numberColumns; iColumn++) { |
---|
1031 | if (!columnLength[iColumn]) |
---|
1032 | continue; |
---|
1033 | if (columnLower[iColumn] < 0.0 || columnUpper[iColumn] > 1.0) |
---|
1034 | good = false; |
---|
1035 | CoinBigIndex j; |
---|
1036 | int nSOS=0; |
---|
1037 | int iSOS=-1; |
---|
1038 | if (!solver->isInteger(iColumn)) |
---|
1039 | good = false; |
---|
1040 | for (j = columnStart[iColumn]; |
---|
1041 | j < columnStart[iColumn] + columnLength[iColumn]; j++) { |
---|
1042 | if (element[j] < 0.0) |
---|
1043 | good = false; |
---|
1044 | int iRow = row[j]; |
---|
1045 | if (rhs[iRow]==-1.0 && sos[iRow] == 'a') { |
---|
1046 | if (element[j] != 1.0) |
---|
1047 | good = false; |
---|
1048 | iSOS=iRow; |
---|
1049 | nSOS++; |
---|
1050 | } |
---|
1051 | } |
---|
1052 | if (nSOS>1) |
---|
1053 | good = false; |
---|
1054 | else if (!nSOS) |
---|
1055 | nonSOS++; |
---|
1056 | sosRow[iColumn] = iSOS; |
---|
1057 | } |
---|
1058 | if (!good) { |
---|
1059 | delete [] sosRow; |
---|
1060 | delete [] rhs; |
---|
1061 | delete [] sos; |
---|
1062 | setWhen(0); // switch off |
---|
1063 | return 0; |
---|
1064 | } |
---|
1065 | } else { |
---|
1066 | abort(); // not allowed yet |
---|
1067 | } |
---|
1068 | const double * solution = solver->getColSolution(); |
---|
1069 | const double * objective = solver->getObjCoefficients(); |
---|
1070 | const double * rowLower = solver->getRowLower(); |
---|
1071 | const double * rowUpper = solver->getRowUpper(); |
---|
1072 | double integerTolerance = model_->getDblParam(CbcModel::CbcIntegerTolerance); |
---|
1073 | double primalTolerance; |
---|
1074 | solver->getDblParam(OsiPrimalTolerance, primalTolerance); |
---|
1075 | |
---|
1076 | numRuns_++; |
---|
1077 | assert (numberRows == matrix_.getNumRows()); |
---|
1078 | // set up linked list for sets |
---|
1079 | int * firstGub = new int [numberRows]; |
---|
1080 | int * nextGub = new int [numberColumns]; |
---|
1081 | int iRow, iColumn; |
---|
1082 | double direction = solver->getObjSense(); |
---|
1083 | double * slackCost = new double [numberRows]; |
---|
1084 | double * modifiedCost = CoinCopyOfArray(objective,numberColumns); |
---|
1085 | for (int iRow = 0;iRow < numberRows; iRow++) { |
---|
1086 | slackCost[iRow]=1.0e30; |
---|
1087 | firstGub[iRow]=-1; |
---|
1088 | } |
---|
1089 | // Take off cost of gub slack |
---|
1090 | for (int iColumn = 0; iColumn < numberColumns; iColumn++) { |
---|
1091 | nextGub[iColumn]=-1; |
---|
1092 | int iRow = sosRow[iColumn]; |
---|
1093 | if (columnLength[iColumn] == 1&&iRow>=0) { |
---|
1094 | // SOS slack |
---|
1095 | double cost = direction*objective[iColumn]; |
---|
1096 | assert (rhs[iRow]<0.0); |
---|
1097 | slackCost[iRow]=CoinMin(slackCost[iRow],cost); |
---|
1098 | } |
---|
1099 | } |
---|
1100 | double offset2 = 0.0; |
---|
1101 | for (int iRow = 0;iRow < numberRows; iRow++) { |
---|
1102 | if (sos[iRow]=='a') { |
---|
1103 | // row is possible |
---|
1104 | sos[iRow]=0; |
---|
1105 | if (rhs[iRow]<0.0) { |
---|
1106 | sos[iRow]=1; |
---|
1107 | rhs[iRow]=1.0; |
---|
1108 | } else if (rhs[iRow] != rowUpper[iRow]) { |
---|
1109 | // G row |
---|
1110 | sos[iRow]=-1; |
---|
1111 | } |
---|
1112 | if( slackCost[iRow] == 1.0e30) { |
---|
1113 | slackCost[iRow]=0.0; |
---|
1114 | } else { |
---|
1115 | offset2 += slackCost[iRow]; |
---|
1116 | sos[iRow] = 2; |
---|
1117 | } |
---|
1118 | } |
---|
1119 | } |
---|
1120 | for (int iColumn = 0; iColumn < numberColumns; iColumn++) { |
---|
1121 | double cost = direction * modifiedCost[iColumn]; |
---|
1122 | CoinBigIndex j; |
---|
1123 | for (j = columnStart[iColumn]; |
---|
1124 | j < columnStart[iColumn] + columnLength[iColumn]; j++) { |
---|
1125 | int iRow = row[j]; |
---|
1126 | if (sos[iRow]>0) { |
---|
1127 | cost -= slackCost[iRow]; |
---|
1128 | if (firstGub[iRow]<0) { |
---|
1129 | firstGub[iRow]=iColumn; |
---|
1130 | } else { |
---|
1131 | int jColumn = firstGub[iRow]; |
---|
1132 | while (nextGub[jColumn]>=0) |
---|
1133 | jColumn=nextGub[jColumn]; |
---|
1134 | nextGub[jColumn]=iColumn; |
---|
1135 | } |
---|
1136 | // Only in one sos |
---|
1137 | break; |
---|
1138 | } |
---|
1139 | } |
---|
1140 | modifiedCost[iColumn] = cost; |
---|
1141 | } |
---|
1142 | delete [] slackCost; |
---|
1143 | double offset; |
---|
1144 | solver->getDblParam(OsiObjOffset, offset); |
---|
1145 | double newSolutionValue = -offset+offset2; |
---|
1146 | int returnCode = 0; |
---|
1147 | |
---|
1148 | |
---|
1149 | // Get solution array for heuristic solution |
---|
1150 | double * newSolution = new double [numberColumns]; |
---|
1151 | double * rowActivity = new double[numberRows]; |
---|
1152 | double * contribution = new double [numberColumns]; |
---|
1153 | int * which = new int [numberColumns]; |
---|
1154 | double * newSolution0 = new double [numberColumns]; |
---|
1155 | if ((algorithm_&(2|4))==0) { |
---|
1156 | // get solution as small as possible |
---|
1157 | for (iColumn = 0; iColumn < numberColumns; iColumn++) |
---|
1158 | newSolution0[iColumn] = columnLower[iColumn]; |
---|
1159 | } else { |
---|
1160 | // Get rounded down solution |
---|
1161 | for (iColumn = 0; iColumn < numberColumns; iColumn++) { |
---|
1162 | double value = solution[iColumn]; |
---|
1163 | // Round down integer |
---|
1164 | if (fabs(floor(value + 0.5) - value) < integerTolerance) { |
---|
1165 | value = floor(CoinMax(value + 1.0e-3, columnLower[iColumn])); |
---|
1166 | } else { |
---|
1167 | value = CoinMax(floor(value), columnLower[iColumn]); |
---|
1168 | } |
---|
1169 | // make sure clean |
---|
1170 | value = CoinMin(value, columnUpper[iColumn]); |
---|
1171 | value = CoinMax(value, columnLower[iColumn]); |
---|
1172 | newSolution0[iColumn] = value; |
---|
1173 | } |
---|
1174 | } |
---|
1175 | double * rowWeight = new double [numberRows]; |
---|
1176 | for (int i=0;i<numberRows;i++) |
---|
1177 | rowWeight[i]=1.0; |
---|
1178 | double costBias = 0.0; |
---|
1179 | int nPass = ((algorithm_&4)!=0) ? 1 : 10; |
---|
1180 | for (int iPass=0;iPass<nPass;iPass++) { |
---|
1181 | newSolutionValue = -offset+offset2; |
---|
1182 | memcpy(newSolution,newSolution0,numberColumns*sizeof(double)); |
---|
1183 | // get row activity |
---|
1184 | memset(rowActivity, 0, numberRows*sizeof(double)); |
---|
1185 | for (iColumn = 0; iColumn < numberColumns; iColumn++) { |
---|
1186 | CoinBigIndex j; |
---|
1187 | double value = newSolution[iColumn]; |
---|
1188 | for (j = columnStart[iColumn]; |
---|
1189 | j < columnStart[iColumn] + columnLength[iColumn]; j++) { |
---|
1190 | int iRow = row[j]; |
---|
1191 | rowActivity[iRow] += value * element[j]; |
---|
1192 | } |
---|
1193 | } |
---|
1194 | if (!rowWeight) { |
---|
1195 | rowWeight = CoinCopyOfArray(rowActivity,numberRows); |
---|
1196 | } |
---|
1197 | for (iColumn = 0; iColumn < numberColumns; iColumn++) { |
---|
1198 | CoinBigIndex j; |
---|
1199 | double value = newSolution[iColumn]; |
---|
1200 | double cost = modifiedCost[iColumn]; |
---|
1201 | double forSort = 1.0e-24; |
---|
1202 | bool hasSlack=false; |
---|
1203 | bool willFit=true; |
---|
1204 | bool gRow=false; |
---|
1205 | newSolutionValue += value * cost; |
---|
1206 | cost += 1.0e-12; |
---|
1207 | for (j = columnStart[iColumn]; |
---|
1208 | j < columnStart[iColumn] + columnLength[iColumn]; j++) { |
---|
1209 | int iRow = row[j]; |
---|
1210 | int type = sos[iRow]; |
---|
1211 | double gap = rhs[iRow] - rowActivity[iRow]+1.0e-8; |
---|
1212 | switch (type) { |
---|
1213 | case -1: |
---|
1214 | // G row |
---|
1215 | gRow = true; |
---|
1216 | #if 0 |
---|
1217 | if (rhs[iRow]>rowWeight[iRow]||(algorithm_&(2|4))!=0) |
---|
1218 | forSort += element[j]; |
---|
1219 | else |
---|
1220 | forSort += 0.1*element[j]; |
---|
1221 | #else |
---|
1222 | forSort += rowWeight[iRow]*element[j]; |
---|
1223 | #endif |
---|
1224 | break; |
---|
1225 | case 0: |
---|
1226 | // L row |
---|
1227 | if (gap<element[j]) { |
---|
1228 | willFit = false; |
---|
1229 | } else { |
---|
1230 | forSort += element[j]; |
---|
1231 | } |
---|
1232 | break; |
---|
1233 | case 1: |
---|
1234 | // SOS without slack |
---|
1235 | if (gap<element[j]) { |
---|
1236 | willFit = false; |
---|
1237 | } |
---|
1238 | break; |
---|
1239 | case 2: |
---|
1240 | // SOS with slack |
---|
1241 | hasSlack = true; |
---|
1242 | if (gap<element[j]) { |
---|
1243 | willFit = false; |
---|
1244 | } |
---|
1245 | break; |
---|
1246 | } |
---|
1247 | } |
---|
1248 | bool isSlack = hasSlack && (columnLength[iColumn]==1); |
---|
1249 | if (forSort<1.0e-24) |
---|
1250 | forSort = 1.0e-12; |
---|
1251 | if ((algorithm_&4)!=0 && forSort > 1.0e-24) |
---|
1252 | forSort=1.0; |
---|
1253 | // Use smallest cost if will fit |
---|
1254 | if (willFit && (hasSlack||gRow) && |
---|
1255 | value == 0.0 && columnUpper[iColumn]) { |
---|
1256 | if (hasSlack && !gRow) { |
---|
1257 | if (cost>1.0e-12) { |
---|
1258 | forSort = 2.0e30; |
---|
1259 | } else if (cost==1.0e-12) { |
---|
1260 | if (!isSlack) |
---|
1261 | forSort = 1.0e29; |
---|
1262 | else |
---|
1263 | forSort = 1.0e28; |
---|
1264 | } else { |
---|
1265 | forSort = cost/forSort; |
---|
1266 | } |
---|
1267 | } else { |
---|
1268 | if (!gRow||true) |
---|
1269 | forSort = (cost+costBias)/forSort; |
---|
1270 | else |
---|
1271 | forSort = 1.0e-12/forSort; |
---|
1272 | } |
---|
1273 | } else { |
---|
1274 | // put at end |
---|
1275 | forSort = 1.0e30; |
---|
1276 | } |
---|
1277 | which[iColumn]=iColumn; |
---|
1278 | contribution[iColumn]= forSort; |
---|
1279 | } |
---|
1280 | CoinSort_2(contribution,contribution+numberColumns,which); |
---|
1281 | // Go through columns |
---|
1282 | int nAdded=0; |
---|
1283 | int nSlacks=0; |
---|
1284 | for (int jColumn = 0; jColumn < numberColumns; jColumn++) { |
---|
1285 | if (contribution[jColumn]>=1.0e30) |
---|
1286 | break; |
---|
1287 | int iColumn = which[jColumn]; |
---|
1288 | double value = newSolution[iColumn]; |
---|
1289 | if (value) |
---|
1290 | continue; |
---|
1291 | bool possible = true; |
---|
1292 | CoinBigIndex j; |
---|
1293 | for (j = columnStart[iColumn]; |
---|
1294 | j < columnStart[iColumn] + columnLength[iColumn]; j++) { |
---|
1295 | int iRow = row[j]; |
---|
1296 | if (sos[iRow]>0&&rowActivity[iRow]) { |
---|
1297 | possible = false; |
---|
1298 | } else { |
---|
1299 | double gap = rhs[iRow] - rowActivity[iRow]+1.0e-8; |
---|
1300 | if (gap<element[j]&&sos[iRow]>=0) |
---|
1301 | possible = false; |
---|
1302 | } |
---|
1303 | } |
---|
1304 | if (possible) { |
---|
1305 | //#define REPORT 1 |
---|
1306 | #ifdef REPORT |
---|
1307 | if ((nAdded%1000)==0) { |
---|
1308 | double gap=0.0; |
---|
1309 | for (int i=0;i<numberRows;i++) { |
---|
1310 | if (rowUpper[i]>1.0e20) |
---|
1311 | gap += CoinMax(rowLower[i]-rowActivity[i],0.0); |
---|
1312 | } |
---|
1313 | if (gap) |
---|
1314 | printf("after %d added gap %g - %d slacks\n", |
---|
1315 | nAdded,gap,nSlacks); |
---|
1316 | } |
---|
1317 | #endif |
---|
1318 | nAdded++; |
---|
1319 | if (columnLength[iColumn]==1) |
---|
1320 | nSlacks++; |
---|
1321 | // Increase chosen column |
---|
1322 | newSolution[iColumn] = 1.0; |
---|
1323 | double cost = modifiedCost[iColumn]; |
---|
1324 | newSolutionValue += cost; |
---|
1325 | for (CoinBigIndex j = columnStart[iColumn]; |
---|
1326 | j < columnStart[iColumn] + columnLength[iColumn]; j++) { |
---|
1327 | int iRow = row[j]; |
---|
1328 | rowActivity[iRow] += element[j]; |
---|
1329 | } |
---|
1330 | } |
---|
1331 | } |
---|
1332 | #ifdef REPORT |
---|
1333 | { |
---|
1334 | double under=0.0; |
---|
1335 | double over=0.0; |
---|
1336 | double gap = 0.0; |
---|
1337 | int nUnder=0; |
---|
1338 | int nOver=0; |
---|
1339 | int nGap=0; |
---|
1340 | for (iRow = 0; iRow < numberRows; iRow++) { |
---|
1341 | if (rowActivity[iRow] < rowLower[iRow] - 10.0*primalTolerance) { |
---|
1342 | double value = rowLower[iRow]-rowActivity[iRow]; |
---|
1343 | #if REPORT>1 |
---|
1344 | printf("below on %d is %g - activity %g lower %g\n", |
---|
1345 | iRow,value,rowActivity[iRow],rowLower[iRow]); |
---|
1346 | #endif |
---|
1347 | under += value; |
---|
1348 | nUnder++; |
---|
1349 | } else if (rowActivity[iRow] > rowUpper[iRow] + 10.0*primalTolerance) { |
---|
1350 | double value = rowActivity[iRow]-rowUpper[iRow]; |
---|
1351 | #if REPORT>1 |
---|
1352 | printf("above on %d is %g - activity %g upper %g\n", |
---|
1353 | iRow,value,rowActivity[iRow],rowUpper[iRow]); |
---|
1354 | #endif |
---|
1355 | over += value; |
---|
1356 | nOver++; |
---|
1357 | } else { |
---|
1358 | double value = rowActivity[iRow]-rowLower[iRow]; |
---|
1359 | if (value && value < 1.0e20) { |
---|
1360 | #if REPORT>1 |
---|
1361 | printf("gap on %d is %g - activity %g lower %g\n", |
---|
1362 | iRow,value,rowActivity[iRow],rowLower[iRow]); |
---|
1363 | #endif |
---|
1364 | gap += value; |
---|
1365 | nGap++; |
---|
1366 | } |
---|
1367 | } |
---|
1368 | } |
---|
1369 | printf("final under %g (%d) - over %g (%d) - free %g (%d) - %d added - solvalue %g\n", |
---|
1370 | under,nUnder,over,nOver,gap,nGap,nAdded,newSolutionValue); |
---|
1371 | } |
---|
1372 | #endif |
---|
1373 | double gap = 0.0; |
---|
1374 | double over = 0.0; |
---|
1375 | int nL=0; |
---|
1376 | int nG=0; |
---|
1377 | int nUnder=0; |
---|
1378 | for (iRow = 0; iRow < numberRows; iRow++) { |
---|
1379 | if (rowLower[iRow]<-1.0e20) |
---|
1380 | nL++; |
---|
1381 | if (rowUpper[iRow]>1.0e20) |
---|
1382 | nG++; |
---|
1383 | if (rowActivity[iRow] < rowLower[iRow] - 10.0*primalTolerance) { |
---|
1384 | gap += rowLower[iRow]-rowActivity[iRow]; |
---|
1385 | nUnder++; |
---|
1386 | rowWeight[iRow] *= 1.1; |
---|
1387 | } else if (rowActivity[iRow] > rowUpper[iRow] + 10.0*primalTolerance) { |
---|
1388 | gap += rowActivity[iRow]-rowUpper[iRow]; |
---|
1389 | } else { |
---|
1390 | over += rowActivity[iRow]-rowLower[iRow]; |
---|
1391 | //rowWeight[iRow] *= 0.9; |
---|
1392 | } |
---|
1393 | } |
---|
1394 | if (nG&&!nL) { |
---|
1395 | // can we fix |
---|
1396 | // get list of columns which can go down without making |
---|
1397 | // things much worse |
---|
1398 | int nPossible=0; |
---|
1399 | int nEasyDown=0; |
---|
1400 | int nSlackDown=0; |
---|
1401 | for (int iColumn=0;iColumn<numberColumns;iColumn++) { |
---|
1402 | if (newSolution[iColumn]&& |
---|
1403 | columnUpper[iColumn]>columnLower[iColumn]) { |
---|
1404 | bool canGoDown=true; |
---|
1405 | bool under = false; |
---|
1406 | int iSos=-1; |
---|
1407 | for (CoinBigIndex j = columnStart[iColumn]; |
---|
1408 | j < columnStart[iColumn] + columnLength[iColumn]; j++) { |
---|
1409 | int iRow = row[j]; |
---|
1410 | if (sos[iRow]<0) { |
---|
1411 | double over = rowActivity[iRow]-rowLower[iRow]; |
---|
1412 | if (over>=0.0&&element[j]>over+1.0e-12) { |
---|
1413 | canGoDown=false; |
---|
1414 | break; |
---|
1415 | } else if (over<0.0) { |
---|
1416 | under = true; |
---|
1417 | } |
---|
1418 | } else { |
---|
1419 | iSos=iRow; |
---|
1420 | } |
---|
1421 | } |
---|
1422 | if (canGoDown) { |
---|
1423 | if (!under) { |
---|
1424 | if (iSos>=0) { |
---|
1425 | // find cheapest |
---|
1426 | double cheapest=modifiedCost[iColumn]; |
---|
1427 | int iCheapest = -1; |
---|
1428 | int jColumn = firstGub[iSos]; |
---|
1429 | assert (jColumn>=0); |
---|
1430 | while (jColumn>=0) { |
---|
1431 | if (modifiedCost[jColumn]<cheapest) { |
---|
1432 | cheapest=modifiedCost[jColumn]; |
---|
1433 | iCheapest=jColumn; |
---|
1434 | } |
---|
1435 | jColumn = nextGub[jColumn]; |
---|
1436 | } |
---|
1437 | if (iCheapest>=0) { |
---|
1438 | // Decrease column |
---|
1439 | newSolution[iColumn] = 0.0; |
---|
1440 | newSolutionValue -= modifiedCost[iColumn]; |
---|
1441 | for (CoinBigIndex j = columnStart[iColumn]; |
---|
1442 | j < columnStart[iColumn] + columnLength[iColumn]; j++) { |
---|
1443 | int iRow = row[j]; |
---|
1444 | rowActivity[iRow] -= element[j]; |
---|
1445 | } |
---|
1446 | // Increase chosen column |
---|
1447 | newSolution[iCheapest] = 1.0; |
---|
1448 | newSolutionValue += modifiedCost[iCheapest]; |
---|
1449 | for (CoinBigIndex j = columnStart[iCheapest]; |
---|
1450 | j < columnStart[iCheapest] + columnLength[iCheapest]; j++) { |
---|
1451 | int iRow = row[j]; |
---|
1452 | rowActivity[iRow] += element[j]; |
---|
1453 | } |
---|
1454 | nEasyDown++; |
---|
1455 | if (columnLength[iColumn]>1) { |
---|
1456 | //printf("%d is easy down\n",iColumn); |
---|
1457 | } else { |
---|
1458 | nSlackDown++; |
---|
1459 | } |
---|
1460 | } |
---|
1461 | } else if (modifiedCost[iColumn]>0.0) { |
---|
1462 | // easy down |
---|
1463 | // Decrease column |
---|
1464 | newSolution[iColumn] = 0.0; |
---|
1465 | newSolutionValue -= modifiedCost[iColumn]; |
---|
1466 | for (CoinBigIndex j = columnStart[iColumn]; |
---|
1467 | j < columnStart[iColumn] + columnLength[iColumn]; j++) { |
---|
1468 | int iRow = row[j]; |
---|
1469 | rowActivity[iRow] -= element[j]; |
---|
1470 | } |
---|
1471 | nEasyDown++; |
---|
1472 | } |
---|
1473 | } else { |
---|
1474 | which[nPossible++]=iColumn; |
---|
1475 | } |
---|
1476 | } |
---|
1477 | } |
---|
1478 | } |
---|
1479 | #ifdef REPORT |
---|
1480 | printf("%d possible down, %d easy down of which %d are slacks\n", |
---|
1481 | nPossible,nEasyDown,nSlackDown); |
---|
1482 | #endif |
---|
1483 | double * needed = new double [numberRows]; |
---|
1484 | for (int i=0;i<numberRows;i++) { |
---|
1485 | double value = rowLower[i] - rowActivity[i]; |
---|
1486 | if (value<1.0e-8) |
---|
1487 | value=0.0; |
---|
1488 | needed[i]=value; |
---|
1489 | } |
---|
1490 | if (gap && /*nUnder==1 &&*/ nonSOS) { |
---|
1491 | double * weight = new double [numberColumns]; |
---|
1492 | int * sort = new int [numberColumns]; |
---|
1493 | // look at ones not in set |
---|
1494 | int nPossible=0; |
---|
1495 | for (int iColumn=0;iColumn<numberColumns;iColumn++) { |
---|
1496 | if (!newSolution[iColumn]&& |
---|
1497 | columnUpper[iColumn]>columnLower[iColumn]) { |
---|
1498 | int iSos=-1; |
---|
1499 | double value=0.0; |
---|
1500 | for (CoinBigIndex j = columnStart[iColumn]; |
---|
1501 | j < columnStart[iColumn] + columnLength[iColumn]; j++) { |
---|
1502 | int iRow = row[j]; |
---|
1503 | if (sos[iRow]<0) { |
---|
1504 | if (needed[iRow]) |
---|
1505 | value += CoinMin(element[j]/needed[iRow],1.0); |
---|
1506 | } else { |
---|
1507 | iSos=iRow; |
---|
1508 | } |
---|
1509 | } |
---|
1510 | if (value && iSos<0) { |
---|
1511 | weight[nPossible]=-value; |
---|
1512 | sort[nPossible++]=iColumn; |
---|
1513 | } |
---|
1514 | } |
---|
1515 | } |
---|
1516 | CoinSort_2(weight,weight+nPossible,sort); |
---|
1517 | for (int i=0;i<nPossible;i++) { |
---|
1518 | int iColumn = sort[i]; |
---|
1519 | double helps=0.0; |
---|
1520 | for (CoinBigIndex j = columnStart[iColumn]; |
---|
1521 | j < columnStart[iColumn] + columnLength[iColumn]; j++) { |
---|
1522 | int iRow = row[j]; |
---|
1523 | if (needed[iRow]) |
---|
1524 | helps += CoinMin(needed[iRow],element[j]); |
---|
1525 | } |
---|
1526 | if (helps) { |
---|
1527 | newSolution[iColumn] = 1.0; |
---|
1528 | newSolutionValue += modifiedCost[iColumn]; |
---|
1529 | for (CoinBigIndex j = columnStart[iColumn]; |
---|
1530 | j < columnStart[iColumn] + columnLength[iColumn]; j++) { |
---|
1531 | int iRow = row[j]; |
---|
1532 | rowActivity[iRow] += element[j]; |
---|
1533 | if (needed[iRow]) { |
---|
1534 | needed[iRow] -= element[j]; |
---|
1535 | if (needed[iRow]<1.0e-8) |
---|
1536 | needed[iRow]=0.0; |
---|
1537 | } |
---|
1538 | } |
---|
1539 | gap -= helps; |
---|
1540 | #ifdef REPORT |
---|
1541 | { |
---|
1542 | double gap2 = 0.0; |
---|
1543 | for (iRow = 0; iRow < numberRows; iRow++) { |
---|
1544 | if (rowActivity[iRow] < rowLower[iRow] - 10.0*primalTolerance) { |
---|
1545 | gap2 += rowLower[iRow]-rowActivity[iRow]; |
---|
1546 | } |
---|
1547 | } |
---|
1548 | printf("estimated gap (nonsos) %g - computed %g\n", |
---|
1549 | gap,gap2); |
---|
1550 | } |
---|
1551 | #endif |
---|
1552 | if (gap<1.0e-12) |
---|
1553 | break; |
---|
1554 | } |
---|
1555 | } |
---|
1556 | delete [] weight; |
---|
1557 | delete [] sort; |
---|
1558 | } |
---|
1559 | if (gap&&nPossible/*&&nUnder==1*/&&true&&model_->bestSolution()) { |
---|
1560 | double * weight = new double [numberColumns]; |
---|
1561 | int * sort = new int [numberColumns]; |
---|
1562 | // look at ones in sets |
---|
1563 | const double * goodSolution = model_->bestSolution(); |
---|
1564 | int nPossible=0; |
---|
1565 | double largestWeight=0.0; |
---|
1566 | for (int iColumn=0;iColumn<numberColumns;iColumn++) { |
---|
1567 | if (!newSolution[iColumn]&&goodSolution[iColumn]&& |
---|
1568 | columnUpper[iColumn]>columnLower[iColumn]) { |
---|
1569 | int iSos=-1; |
---|
1570 | double value=0.0; |
---|
1571 | for (CoinBigIndex j = columnStart[iColumn]; |
---|
1572 | j < columnStart[iColumn] + columnLength[iColumn]; j++) { |
---|
1573 | int iRow = row[j]; |
---|
1574 | if (sos[iRow]<0) { |
---|
1575 | if (needed[iRow]) |
---|
1576 | value += CoinMin(element[j]/needed[iRow],1.0); |
---|
1577 | } else { |
---|
1578 | iSos=iRow; |
---|
1579 | } |
---|
1580 | } |
---|
1581 | if (value&&iSos>=0) { |
---|
1582 | // see if value bigger than current |
---|
1583 | int jColumn = firstGub[iSos]; |
---|
1584 | assert (jColumn>=0); |
---|
1585 | while (jColumn>=0) { |
---|
1586 | if (newSolution[jColumn]) |
---|
1587 | break; |
---|
1588 | jColumn = nextGub[jColumn]; |
---|
1589 | } |
---|
1590 | assert (jColumn>=0); |
---|
1591 | double value2=0.0; |
---|
1592 | for (CoinBigIndex j = columnStart[jColumn]; |
---|
1593 | j < columnStart[jColumn] + columnLength[jColumn]; j++) { |
---|
1594 | int iRow = row[j]; |
---|
1595 | if (needed[iRow]) |
---|
1596 | value2 += CoinMin(element[j]/needed[iRow],1.0); |
---|
1597 | } |
---|
1598 | if (value>value2) { |
---|
1599 | weight[nPossible]=-(value-value2); |
---|
1600 | largestWeight = CoinMax(largestWeight,(value-value2)); |
---|
1601 | sort[nPossible++]=iColumn; |
---|
1602 | } |
---|
1603 | } |
---|
1604 | } |
---|
1605 | } |
---|
1606 | if (nPossible) { |
---|
1607 | double * temp = new double [numberRows]; |
---|
1608 | int * which2 = new int [numberRows]; |
---|
1609 | memset(temp,0,numberRows*sizeof(double)); |
---|
1610 | // modify so ones just more than gap best |
---|
1611 | if (largestWeight>gap&&nUnder==1) { |
---|
1612 | double offset = 4*largestWeight; |
---|
1613 | for (int i=0;i<nPossible;i++) { |
---|
1614 | double value = -weight[i]; |
---|
1615 | if (value>gap-1.0e-12) |
---|
1616 | weight[i] = -(offset-(value-gap)); |
---|
1617 | } |
---|
1618 | } |
---|
1619 | CoinSort_2(weight,weight+nPossible,sort); |
---|
1620 | for (int i=0;i<nPossible;i++) { |
---|
1621 | int iColumn = sort[i]; |
---|
1622 | int n=0; |
---|
1623 | // find jColumn |
---|
1624 | int iSos=-1; |
---|
1625 | for (CoinBigIndex j = columnStart[iColumn]; |
---|
1626 | j < columnStart[iColumn] + columnLength[iColumn]; j++) { |
---|
1627 | int iRow = row[j]; |
---|
1628 | temp[iRow]=element[j]; |
---|
1629 | which2[n++]=iRow; |
---|
1630 | if (sos[iRow]>=0) { |
---|
1631 | iSos=iRow; |
---|
1632 | } |
---|
1633 | } |
---|
1634 | int jColumn = firstGub[iSos]; |
---|
1635 | assert (jColumn>=0); |
---|
1636 | while (jColumn>=0) { |
---|
1637 | if (newSolution[jColumn]) |
---|
1638 | break; |
---|
1639 | jColumn = nextGub[jColumn]; |
---|
1640 | } |
---|
1641 | assert (jColumn>=0); |
---|
1642 | for (CoinBigIndex j = columnStart[jColumn]; |
---|
1643 | j < columnStart[jColumn] + columnLength[jColumn]; j++) { |
---|
1644 | int iRow = row[j]; |
---|
1645 | if (!temp[iRow]) |
---|
1646 | which2[n++]=iRow; |
---|
1647 | temp[iRow] -= element[j]; |
---|
1648 | } |
---|
1649 | double helps = 0.0; |
---|
1650 | for (int i=0;i<n;i++) { |
---|
1651 | int iRow = which2[i]; |
---|
1652 | double newValue = rowActivity[iRow]+temp[iRow]; |
---|
1653 | if (temp[iRow]>1.0e-8) { |
---|
1654 | if (rowActivity[iRow]<rowLower[iRow]-1.0e-8) { |
---|
1655 | helps += CoinMin(temp[iRow], |
---|
1656 | rowLower[iRow]-rowActivity[iRow]); |
---|
1657 | } |
---|
1658 | } else if (temp[iRow]<-1.0e-8) { |
---|
1659 | if (newValue<rowLower[iRow]-1.0e-12) { |
---|
1660 | helps -= CoinMin(-temp[iRow], |
---|
1661 | 1.0*(rowLower[iRow]-newValue)); |
---|
1662 | } |
---|
1663 | } |
---|
1664 | } |
---|
1665 | if (helps>0.0) { |
---|
1666 | newSolution[iColumn]=1.0; |
---|
1667 | newSolution[jColumn]=0.0; |
---|
1668 | newSolutionValue += modifiedCost[iColumn]-modifiedCost[jColumn]; |
---|
1669 | for (int i=0;i<n;i++) { |
---|
1670 | int iRow = which2[i]; |
---|
1671 | double newValue = rowActivity[iRow]+temp[iRow]; |
---|
1672 | rowActivity[iRow] = newValue; |
---|
1673 | temp[iRow]=0.0; |
---|
1674 | } |
---|
1675 | gap -= helps; |
---|
1676 | #ifdef REPORT |
---|
1677 | { |
---|
1678 | double gap2 = 0.0; |
---|
1679 | for (iRow = 0; iRow < numberRows; iRow++) { |
---|
1680 | if (rowActivity[iRow] < rowLower[iRow] - 10.0*primalTolerance) { |
---|
1681 | gap2 += rowLower[iRow]-rowActivity[iRow]; |
---|
1682 | } |
---|
1683 | } |
---|
1684 | printf("estimated gap %g - computed %g\n", |
---|
1685 | gap,gap2); |
---|
1686 | } |
---|
1687 | #endif |
---|
1688 | if (gap<1.0e-8) |
---|
1689 | break; |
---|
1690 | } else { |
---|
1691 | for (int i=0;i<n;i++) |
---|
1692 | temp[which2[i]]=0.0; |
---|
1693 | } |
---|
1694 | } |
---|
1695 | delete [] which2; |
---|
1696 | delete [] temp; |
---|
1697 | } |
---|
1698 | delete [] weight; |
---|
1699 | delete [] sort; |
---|
1700 | } |
---|
1701 | delete [] needed; |
---|
1702 | } |
---|
1703 | #ifdef REPORT |
---|
1704 | { |
---|
1705 | double gap=0.0; |
---|
1706 | double over = 0.0; |
---|
1707 | for (iRow = 0; iRow < numberRows; iRow++) { |
---|
1708 | if (rowActivity[iRow] < rowLower[iRow] - 10.0*primalTolerance) { |
---|
1709 | double value = rowLower[iRow]-rowActivity[iRow]; |
---|
1710 | #if REPORT>1 |
---|
1711 | printf("below on %d is %g - activity %g lower %g\n", |
---|
1712 | iRow,value,rowActivity[iRow],rowLower[iRow]); |
---|
1713 | #endif |
---|
1714 | gap += value; |
---|
1715 | } else if (rowActivity[iRow] > rowUpper[iRow] + 10.0*primalTolerance) { |
---|
1716 | double value = rowActivity[iRow]-rowUpper[iRow]; |
---|
1717 | #if REPORT>1 |
---|
1718 | printf("above on %d is %g - activity %g upper %g\n", |
---|
1719 | iRow,value,rowActivity[iRow],rowUpper[iRow]); |
---|
1720 | #endif |
---|
1721 | gap += value; |
---|
1722 | } else { |
---|
1723 | double value = rowActivity[iRow]-rowLower[iRow]; |
---|
1724 | if (value) { |
---|
1725 | #if REPORT>1 |
---|
1726 | printf("over on %d is %g - activity %g lower %g\n", |
---|
1727 | iRow,value,rowActivity[iRow],rowLower[iRow]); |
---|
1728 | #endif |
---|
1729 | over += value; |
---|
1730 | } |
---|
1731 | } |
---|
1732 | } |
---|
1733 | printf("modified final gap %g - over %g - %d added - solvalue %g\n", |
---|
1734 | gap,over,nAdded,newSolutionValue); |
---|
1735 | } |
---|
1736 | #endif |
---|
1737 | if (!gap) { |
---|
1738 | break; |
---|
1739 | } else { |
---|
1740 | if (iPass==0) { |
---|
1741 | costBias = 10.0*newSolutionValue/static_cast<double>(nAdded); |
---|
1742 | } else { |
---|
1743 | costBias *= 10.0; |
---|
1744 | } |
---|
1745 | } |
---|
1746 | } |
---|
1747 | delete [] newSolution0; |
---|
1748 | delete [] rowWeight; |
---|
1749 | delete [] sos; |
---|
1750 | delete [] firstGub; |
---|
1751 | delete [] nextGub; |
---|
1752 | if (newSolutionValue < solutionValue) { |
---|
1753 | // check feasible |
---|
1754 | memset(rowActivity, 0, numberRows*sizeof(double)); |
---|
1755 | for (iColumn = 0; iColumn < numberColumns; iColumn++) { |
---|
1756 | CoinBigIndex j; |
---|
1757 | double value = newSolution[iColumn]; |
---|
1758 | if (value) { |
---|
1759 | for (j = columnStart[iColumn]; |
---|
1760 | j < columnStart[iColumn] + columnLength[iColumn]; j++) { |
---|
1761 | int iRow = row[j]; |
---|
1762 | rowActivity[iRow] += value * element[j]; |
---|
1763 | } |
---|
1764 | } |
---|
1765 | } |
---|
1766 | // check was approximately feasible |
---|
1767 | bool feasible = true; |
---|
1768 | for (iRow = 0; iRow < numberRows; iRow++) { |
---|
1769 | if (rowActivity[iRow] < rowLower[iRow]) { |
---|
1770 | if (rowActivity[iRow] < rowLower[iRow] - 10.0*primalTolerance) |
---|
1771 | feasible = false; |
---|
1772 | } else if (rowActivity[iRow] > rowUpper[iRow]) { |
---|
1773 | if (rowActivity[iRow] > rowUpper[iRow] + 10.0*primalTolerance) |
---|
1774 | feasible = false; |
---|
1775 | } |
---|
1776 | } |
---|
1777 | if (feasible) { |
---|
1778 | // new solution |
---|
1779 | memcpy(betterSolution, newSolution, numberColumns*sizeof(double)); |
---|
1780 | solutionValue = newSolutionValue; |
---|
1781 | //printf("** Solution of %g found by rounding\n",newSolutionValue); |
---|
1782 | returnCode = 1; |
---|
1783 | } else { |
---|
1784 | // Can easily happen |
---|
1785 | //printf("Debug CbcHeuristicGreedySOS giving bad solution\n"); |
---|
1786 | } |
---|
1787 | } |
---|
1788 | delete [] sosRow; |
---|
1789 | delete [] newSolution; |
---|
1790 | delete [] rowActivity; |
---|
1791 | delete [] modifiedCost; |
---|
1792 | delete [] contribution; |
---|
1793 | delete [] which; |
---|
1794 | delete [] rhs; |
---|
1795 | return returnCode; |
---|
1796 | } |
---|
1797 | // update model |
---|
1798 | void CbcHeuristicGreedySOS::setModel(CbcModel * model) |
---|
1799 | { |
---|
1800 | delete [] originalRhs_; |
---|
1801 | gutsOfConstructor(model); |
---|
1802 | validate(); |
---|
1803 | } |
---|
1804 | // Resets stuff if model changes |
---|
1805 | void |
---|
1806 | CbcHeuristicGreedySOS::resetModel(CbcModel * model) |
---|
1807 | { |
---|
1808 | delete [] originalRhs_; |
---|
1809 | gutsOfConstructor(model); |
---|
1810 | } |
---|
1811 | // Validate model i.e. sets when_ to 0 if necessary (may be NULL) |
---|
1812 | void |
---|
1813 | CbcHeuristicGreedySOS::validate() |
---|
1814 | { |
---|
1815 | if (model_ && when() < 10) { |
---|
1816 | if (model_->numberIntegers() != |
---|
1817 | model_->numberObjects() && (model_->numberObjects() || |
---|
1818 | (model_->specialOptions()&1024) == 0)) { |
---|
1819 | int numberOdd = 0; |
---|
1820 | for (int i = 0; i < model_->numberObjects(); i++) { |
---|
1821 | if (!model_->object(i)->canDoHeuristics()) |
---|
1822 | numberOdd++; |
---|
1823 | } |
---|
1824 | if (numberOdd) |
---|
1825 | setWhen(0); |
---|
1826 | } |
---|
1827 | // Only works if coefficients positive and all rows L/G or SOS |
---|
1828 | OsiSolverInterface * solver = model_->solver(); |
---|
1829 | const double * columnUpper = solver->getColUpper(); |
---|
1830 | const double * columnLower = solver->getColLower(); |
---|
1831 | const double * rowLower = solver->getRowLower(); |
---|
1832 | const double * rowUpper = solver->getRowUpper(); |
---|
1833 | |
---|
1834 | int numberRows = solver->getNumRows(); |
---|
1835 | // Column copy |
---|
1836 | const double * element = matrix_.getElements(); |
---|
1837 | const int * row = matrix_.getIndices(); |
---|
1838 | const CoinBigIndex * columnStart = matrix_.getVectorStarts(); |
---|
1839 | const int * columnLength = matrix_.getVectorLengths(); |
---|
1840 | bool good = true; |
---|
1841 | assert (originalRhs_); |
---|
1842 | for (int iRow = 0; iRow < numberRows; iRow++) { |
---|
1843 | if (rowLower[iRow] == 1.0 && rowUpper[iRow] == 1.0) { |
---|
1844 | // SOS |
---|
1845 | originalRhs_[iRow]=-1.0; |
---|
1846 | } else if (rowLower[iRow] > 0.0 && rowUpper[iRow] < 1.0e10) { |
---|
1847 | good = false; |
---|
1848 | } else if (rowUpper[iRow] < 0.0) { |
---|
1849 | good = false; |
---|
1850 | } else if (rowUpper[iRow] < 1.0e10) { |
---|
1851 | originalRhs_[iRow]=rowUpper[iRow]; |
---|
1852 | } else { |
---|
1853 | originalRhs_[iRow]=rowLower[iRow]; |
---|
1854 | } |
---|
1855 | } |
---|
1856 | int numberColumns = solver->getNumCols(); |
---|
1857 | for (int iColumn = 0; iColumn < numberColumns; iColumn++) { |
---|
1858 | if (!columnLength[iColumn]) |
---|
1859 | continue; |
---|
1860 | if (columnLower[iColumn] < 0.0 || columnUpper[iColumn] > 1.0) |
---|
1861 | good = false; |
---|
1862 | CoinBigIndex j; |
---|
1863 | int nSOS=0; |
---|
1864 | if (!solver->isInteger(iColumn)) |
---|
1865 | good = false; |
---|
1866 | for (j = columnStart[iColumn]; |
---|
1867 | j < columnStart[iColumn] + columnLength[iColumn]; j++) { |
---|
1868 | if (element[j] < 0.0) |
---|
1869 | good = false; |
---|
1870 | int iRow = row[j]; |
---|
1871 | if (originalRhs_[iRow]==-1.0) { |
---|
1872 | if (element[j] != 1.0) |
---|
1873 | good = false; |
---|
1874 | nSOS++; |
---|
1875 | } |
---|
1876 | } |
---|
1877 | if (nSOS > 1) |
---|
1878 | good = false; |
---|
1879 | } |
---|
1880 | if (!good) |
---|
1881 | setWhen(0); // switch off |
---|
1882 | } |
---|
1883 | } |
---|