#
source:
pico/trunk/src/milp/pico/feasibilityPump.cpp
@
5603

Revision 5603, 31.6 KB checked in by jwatson, 6 years ago (diff) |
---|

Line | |
---|---|

1 | /* _________________________________________________________________________ |

2 | * |

3 | * PICO: A C++ library of scalable branch-and-bound and related methods. |

4 | * Copyright (c) 2007, Sandia National Laboratories. |

5 | * This software is distributed under the GNU Lesser General Public License. |

6 | * For more information, see the README file in the top PICO directory. |

7 | * _________________________________________________________________________ |

8 | */ |

9 | |

10 | |

11 | #include <acro_config.h> |

12 | #include <pebbl/gRandom.h> |

13 | #include <pico/milp.h> |

14 | #include <pico/milpNode.h> |

15 | |

16 | #include <pico/feasibilityPump.h> |

17 | |

18 | #include <CoinPackedVector.hpp> |

19 | #include <CoinWarmStart.hpp> |

20 | |

21 | #include <iostream> |

22 | #include <iomanip> |

23 | #include <limits> |

24 | |

25 | #include <math.h> |

26 | #include <assert.h> |

27 | |

28 | // comparison function used in std::sort call found below. |

29 | // sorts in decreasing order. |

30 | static bool compareDists(const std::pair<double,int> & a, |

31 | const std::pair<double,int> & b) |

32 | { |

33 | return (a.first > b.first); |

34 | } |

35 | |

36 | // a pretty-print / debug utility. |

37 | template <typename X> |

38 | void prettyPrint(const char * prefix, |

39 | const X * vector_data, |

40 | unsigned int vector_length, |

41 | std::ostream &ostr, |

42 | int width = 4) |

43 | { |

44 | ostr << prefix << "=" << std::endl; |

45 | ostr << "["; |

46 | for (unsigned int i = 0; i < vector_length; ++i) |

47 | { |

48 | ostr.precision(2); |

49 | ostr << std::setw(width) << vector_data[i]; |

50 | if (i != vector_length - 1) |

51 | { |

52 | ostr << " "; |

53 | } |

54 | } |

55 | ostr << "]" << std::endl; |

56 | } |

57 | |

58 | namespace pico |

59 | { |

60 | |

61 | feasibilityPump::feasibilityPump(void) |

62 | { |

63 | name = "Feasibility pump"; |

64 | |

65 | FPDebug = 0; |

66 | ParameterSet::create_categorized_parameter("FPDebug",FPDebug, |

67 | "<int>","0", |

68 | "Specific debug level for feasibility pump", |

69 | "Debugging", |

70 | utilib::ParameterPositive<int>()); |

71 | |

72 | FPMetaIterations = 1; |

73 | ParameterSet::create_categorized_parameter("FPMetaIterations",FPMetaIterations, |

74 | "<int>","1", |

75 | "Feasibility pump meta-iteration limit", |

76 | "Feasibility Pump Heuristic", |

77 | utilib::ParameterPositive<int>()); |

78 | |

79 | FPRootIterations = 2000; |

80 | ParameterSet::create_categorized_parameter("FPRootIterations",FPRootIterations, |

81 | "<int>","2000", |

82 | "Limit on the cumulative feasibility pump iterations at the root node.", |

83 | "Feasibility Pump Heuristic", |

84 | utilib::ParameterPositive<int>()); |

85 | |

86 | FPInteriorIterations = 50; |

87 | ParameterSet::create_categorized_parameter("FPInteriorIterations",FPInteriorIterations, |

88 | "<int>","50", |

89 | "Limit on the cumulative feasibility pump iterations at interior nodes.", |

90 | "Feasibility Pump Heuristic", |

91 | utilib::ParameterPositive<int>()); |

92 | |

93 | FPRunTime = 300; |

94 | ParameterSet::create_categorized_parameter("FPRunTime",FPRunTime, |

95 | "<double>","300", |

96 | "Limit on the cumulative CPU seconds spent in the feasibility pump heuristic", |

97 | "Feasibility Pump Heuristic", |

98 | utilib::ParameterPositive<double>()); |

99 | |

100 | FPRandomInitialSolution = false; |

101 | ParameterSet::create_categorized_parameter("FPRandomInitialSolution",FPRandomInitialSolution, |

102 | "<bool>","false", |

103 | "Starts feasibilty pump from a random solution", |

104 | "Feasibility Pump Heuristic"); |

105 | |

106 | FPImproveIncumbent = true; |

107 | ParameterSet::create_categorized_parameter("FPImproveIncumbent",FPImproveIncumbent, |

108 | "<bool>","true", |

109 | "Allows feasibility pump to continue inner-loop iterations in attempt to improve an initial incumbent", |

110 | "Feasibility Pump Heuristic"); |

111 | |

112 | FPImprovementAlpha = 0.50; |

113 | ParameterSet::create_categorized_parameter("FPImprovementAlpha",FPImprovementAlpha, |

114 | "<double>","0.50", |

115 | "When improving incumbents, the fraction of the absolute gap used to form the new objective cut", |

116 | "Feasibility Pump Heuristic"); |

117 | |

118 | FPMeanFlipCount = 20; |

119 | ParameterSet::create_categorized_parameter("FPMeanFlipCount",FPMeanFlipCount, |

120 | "<int>","20", |

121 | "Mean number of variables 'flipped' per feasibility pump iteration, if rounded solution remains unchanged", |

122 | "Feasibility Pump Heuristic", |

123 | utilib::ParameterPositive<int>()); |

124 | |

125 | FPRoundingType = 0; |

126 | ParameterSet::create_categorized_parameter("FPRoundingType",FPRoundingType, |

127 | "<int>","0", |

128 | "The type of rounding used by FP to form a new x~ target vector from the current LP solution x*. 0 indicates standard deterministic. 1 indicates randomized rounding around the interval [0.33, 0.66]", |

129 | "Feasibility Pump Heuristic", |

130 | utilib::ParameterPositive<int>()); |

131 | |

132 | FPPerturbationInterval = 100; |

133 | ParameterSet::create_categorized_parameter("FPPerturbationInterval",FPPerturbationInterval, |

134 | "<int>","100", |

135 | "The number of FP iterations between random perturbations of the current x~ solution", |

136 | "Feasibility Pump Heuristic", |

137 | utilib::ParameterPositive<int>()); |

138 | |

139 | } |

140 | |

141 | feasibilityPump::~feasibilityPump(void) |

142 | { |

143 | // no dynamic memory to worry about. |

144 | } |

145 | |

146 | void feasibilityPump::reset(MILP* mGlobal_, lpShareObj* lpShare_) |

147 | { |

148 | lpBasedHeuristic::reset(mGlobal_, lpShare_); |

149 | |

150 | // Serial parameter settings that are hierarchical/not automatic |

151 | |

152 | if (parameter_initialized("FPDebug")) |

153 | { |

154 | setDebug(FPDebug); |

155 | } |

156 | |

157 | // FP needs the objective cut when dealing with |

158 | // non-root-node solves. |

159 | useObjCut = true; |

160 | } |

161 | |

162 | void feasibilityPump::setupLP(MILPNode * sp) |

163 | { |

164 | // we currently execute FP whenever the heuristic management framework tells us to... |

165 | currentDepth = sp->depth; |

166 | executeFP = true; |

167 | |

168 | DEBUGPR(10, ucout << "Setting up feasibility pump to execute at depth=" << sp->depth << std::endl << Flush); |

169 | |

170 | lpBasedHeuristic::setupLP(sp); |

171 | } |

172 | |

173 | MILPSolution * feasibilityPump::runHeuristic(void) |

174 | { |

175 | // execute FP according to whatever rule is encoded in the setupLP method. |

176 | if (executeFP == false) |

177 | { |

178 | return 0; |

179 | } |

180 | |

181 | // stops the LP from thinking that it's just solving the relaxation - |

182 | // we need to know which variables are integer. |

183 | // IMPT: have to un-do on the way out. |

184 | lp->setCutting(); |

185 | |

186 | // verify that the rounding type is set to a sane value. |

187 | if (FPRoundingType >= 2) |

188 | { |

189 | DEBUGPR(10, ucout << "***WARNING: FPRoundingType parameter has an unknown value=" << FPRoundingType << std::endl << Flush); |

190 | return 0; |

191 | } |

192 | |

193 | DEBUGPR(100, ucout << "-----------------------------------------------" << std::endl << Flush); |

194 | DEBUGPR(10, ucout << "Executing feasibility pump" << std::endl << Flush); |

195 | |

196 | double start_time = CPUSeconds(); |

197 | |

198 | int cumulative_iterations(0); |

199 | int allowable_iterations; |

200 | if (currentDepth == 1) // depth is 1-based in PICO/PEBBL. |

201 | { |

202 | allowable_iterations = FPRootIterations; |

203 | } |

204 | else |

205 | { |

206 | allowable_iterations = FPInteriorIterations; |

207 | } |

208 | |

209 | int num_vars(lp->getNumCols()); |

210 | |

211 | // verify that any integer variables are in fact binary - the current |

212 | // code doesn't deal with the general integer case, as that's more |

213 | // than a bit of pain. |

214 | for (int i = 0; i < num_vars; ++i) |

215 | { |

216 | if ((lp->isContinuous(i) == false) && |

217 | (lp->isBinary(i) == false)) |

218 | { |

219 | DEBUGPR(10, ucout << "Not all integer variables are binary - " |

220 | "feasibility pump is currently only coded for binary MIPs" |

221 | << std::endl << Flush); |

222 | return 0; |

223 | } |

224 | } |

225 | |

226 | DEBUGPR(50, ucout << "Limit on number of feasibility pump meta iterations=" |

227 | << FPMetaIterations << std::endl << Flush); |

228 | DEBUGPR(50, ucout << "Limit on cumulative number of feasibility pump " |

229 | "iterations=" << allowable_iterations << std::endl << Flush); |

230 | DEBUGPR(50, ucout << "Limit (in seconds) on the cumulative feasibility pump " |

231 | " cpu time=" << FPRunTime << std::endl << Flush); |

232 | |

233 | // statistics for the best-so-far incumbent. |

234 | double best_objective; |

235 | DoubleVector best_solution(num_vars); |

236 | if (lp->getObjSense() == pebbl::minimize) |

237 | { |

238 | best_objective = std::numeric_limits<double>::max(); |

239 | } |

240 | else |

241 | { |

242 | best_objective = std::numeric_limits<double>::min(); |

243 | } |

244 | |

245 | int executed_iterations(0); |

246 | |

247 | for (int meta_iteration = 1; |

248 | meta_iteration <= FPMetaIterations; |

249 | ++meta_iteration) |

250 | { |

251 | DEBUGPR(10, ucout << "Executing meta-level iteration=" |

252 | << meta_iteration << std::endl << Flush); |

253 | |

254 | ++executed_iterations; |

255 | |

256 | // clone the LP interface to solve the feasibility pump sub-problems. |

257 | // NOTE: if you do this in the outer loop, i.e., create it once, |

258 | // then re-solves with the LP* relaxation initial solution |

259 | // immediately returned with a feasible (high-quality) solution. |

260 | // this wasn't happening for random initial solutions. not |

261 | // resolved yet as to why this happens... |

262 | // NOTE: I moved the clone/resolve to this point, because the |

263 | // lp object has no objective/solution - for reasons we |

264 | // don't understand at all. |

265 | OsiSolverInterface * fp_solver_interface = lp->OsiOnlyCopy(); |

266 | // Cbc does this - isn't clear to me why, but... |

267 | fp_solver_interface->resolve(); |

268 | |

269 | // depending on the instance and the time allocated to FP, the |

270 | // resolve may take more time that is allocated to FP. in this |

271 | // case, break out of the meta-iteration loop and terminate. |

272 | if ((CPUSeconds() - start_time) >= FPRunTime) |

273 | { |

274 | DEBUGPR(10, ucout << "Time limit exceeded during resolve() call to cloned LP solver for FP" << std::endl << Flush); |

275 | break; |

276 | } |

277 | |

278 | // form an integer (but likely infeasible) starting solution. |

279 | double * initial_solution = new double[num_vars]; |

280 | double initial_solution_lb = fp_solver_interface->getInfinity() * fp_solver_interface->getObjSense(); |

281 | |

282 | if (FPRandomInitialSolution == true) |

283 | { |

284 | DEBUGPR(10, ucout << "Starting search from random solution" |

285 | << std::endl << Flush); |

286 | |

287 | RNG * rng = pebbl::gRandomRNG(); |

288 | |

289 | for (int i = 0; i < num_vars; ++i) |

290 | { |

291 | if (lp->isContinuous(i) == true) |

292 | { |

293 | initial_solution[i] = 0.0; // isn't manipulated by FP anyway |

294 | } |

295 | else |

296 | { |

297 | DUniform<int> sampler(rng, |

298 | int(rint(lp->getColLower(i))), |

299 | int(rint(lp->getColUpper(i)))); |

300 | initial_solution[i] = sampler(); |

301 | } |

302 | } |

303 | } |

304 | else |

305 | { |

306 | DEBUGPR(10, ucout << "Starting search from LP relaxation" << std::endl << Flush); |

307 | |

308 | const double * x = fp_solver_interface->getColSolution(); |

309 | |

310 | DEBUGPR(100, prettyPrint("LP relaxation: ", x, num_vars, ucout)); |

311 | DEBUGPR(100, ucout << "Objective=" << fp_solver_interface->getObjValue() << std::endl;); |

312 | |

313 | initial_solution_lb = fp_solver_interface->getObjValue(); |

314 | |

315 | // the first step in the feasibility pump is to solve the |

316 | // LP relaxation, in order to get the initial x^{*} vector. |

317 | // fortunately, we already have that as a side-effect of |

318 | // MILPNode. actually, any solution is OK - the idea behind |

319 | // using the LP relaxation is to have a prayer at finding |

320 | // a low-cost feasible solution. |

321 | |

322 | for (int i = 0; i < num_vars; ++i) |

323 | { |

324 | // if a variable is continuous, leave it alone (it won't |

325 | // be examined anyway). otherwise, round it to the |

326 | // nearest integer value. NOTE: using lp instead of |

327 | // fp_solver_interface because there aren't any integers |

328 | // in the copy for some reason... |

329 | if (lp->isContinuous(i) == true) |

330 | { |

331 | initial_solution[i] = x[i]; // isn't used anyway |

332 | } |

333 | else |

334 | { |

335 | initial_solution[i] = rint(x[i]); |

336 | } |

337 | } |

338 | } |

339 | |

340 | // the feasibility pump objective is always to minimize distance. |

341 | fp_solver_interface->setObjSense(pebbl::minimize); |

342 | |

343 | // add a row to the cloned solver interface to represent |

344 | // the objective cut. at the root node, this starts out at |

345 | // infinity, and is incrementally modified as better and |

346 | // better solutions are found - assuming incumbent improvement |

347 | // is enabled. |

348 | // NOTE: does milpNode or the base heuristic already add the objective cut? |

349 | const double * obj_coeffs = fp_solver_interface->getObjCoefficients(); |

350 | CoinPackedVector objective_cut(num_vars, obj_coeffs); |

351 | fp_solver_interface->addRow(objective_cut, -(fp_solver_interface->getInfinity()), fp_solver_interface->getInfinity()); |

352 | // the new row is added as the last index (num_rows - 1). |

353 | int num_rows(fp_solver_interface->getNumRows()); |

354 | |

355 | // call the core feasibility pump routine. |

356 | double * resulting_solution = new double[num_vars]; |

357 | int iterations_this_trial(0); |

358 | double runtime_this_trial(0.0); |

359 | bool result = feasibilityPumpCore(fp_solver_interface, |

360 | initial_solution, |

361 | initial_solution_lb, |

362 | allowable_iterations |

363 | - cumulative_iterations, |

364 | FPRunTime |

365 | - (CPUSeconds() - start_time), |

366 | iterations_this_trial, |

367 | runtime_this_trial, |

368 | resulting_solution); |

369 | |

370 | cumulative_iterations += iterations_this_trial; |

371 | |

372 | if (result == true) |

373 | { |

374 | DEBUGPR(50, ucout << "Updating incumbent with feasibility pump" |

375 | " solution" << std::endl << Flush); |

376 | |

377 | // translate solution into utilib world. |

378 | DoubleVector heuristic_solution(num_vars); |

379 | for (int i = 0; i < num_vars; ++i) |

380 | { |

381 | heuristic_solution[i] = resulting_solution[i]; |

382 | } |

383 | |

384 | double this_objective = lp->evalVector(heuristic_solution); |

385 | |

386 | if (((lp->getObjSense() == pebbl::minimize) && |

387 | (this_objective < best_objective)) || |

388 | ((lp->getObjSense() == pebbl::maximize) && |

389 | (this_objective > best_objective))) |

390 | { |

391 | best_objective = this_objective; |

392 | best_solution = heuristic_solution; |

393 | } |

394 | } |

395 | else |

396 | { |

397 | break; |

398 | } |

399 | |

400 | // clean up. |

401 | delete [] initial_solution; |

402 | delete [] resulting_solution; |

403 | delete fp_solver_interface; |

404 | } |

405 | |

406 | MILPSolution * result = 0; |

407 | |

408 | if (((lp->getObjSense() == pebbl::minimize) && |

409 | (best_objective != std::numeric_limits<double>::max())) || |

410 | ((lp->getObjSense() == pebbl::maximize) && |

411 | (best_objective != std::numeric_limits<double>::min()))) |

412 | { |

413 | DEBUGPR(10, ucout << "Objective value of best incumbent=" |

414 | << best_objective << ", "); |

415 | DEBUGPR(10, ucout << "Meta-iterations=" << executed_iterations << ", "); |

416 | DEBUGPR(10, ucout << "Cumulative iterations=" << cumulative_iterations |

417 | << ", "); |

418 | DEBUGPR(10, ucout << "Cumulative CPU time=" << CPUSeconds() - start_time |

419 | << " seconds" << std::endl << Flush); |

420 | |

421 | // NOTE: the heuristic solution is unlikely to be exactly integral |

422 | // in the 'integer' variables. may cause problems when in |

423 | // acro validation mode. |

424 | |

425 | result = new MILPSolution(best_solution, mGlobal, lp->evalVector(best_solution)); |

426 | } |

427 | else |

428 | { |

429 | DEBUGPR(10, ucout << "No feasible incumbent was found" << std::endl << Flush); |

430 | } |

431 | |

432 | // revert the LP back to where it was ("solve" mode). |

433 | lp->setSolving(); |

434 | |

435 | DEBUGPR(100, ucout << "-----------------------------------------------" << std::endl << Flush); |

436 | |

437 | |

438 | return result; |

439 | } |

440 | |

441 | bool feasibilityPump::deterministicRounding(unsigned int & num_rounded_up, |

442 | unsigned int & num_rounded_down, |

443 | int num_vars, |

444 | const double * new_solution, |

445 | double * x_tilde)const |

446 | { |

447 | bool any_inverted = false; |

448 | num_rounded_down = 0; |

449 | num_rounded_up = 0; |

450 | for (int j = 0; j < num_vars; ++j) |

451 | { |

452 | if (lp->isContinuous(j) == false) |

453 | { |

454 | if (bGlobal->isInteger(new_solution[j]) == false) |

455 | { |

456 | double rounded_value = rint(new_solution[j]); |

457 | if (fabs(rounded_value - x_tilde[j]) > 1e-5) // enough of a threshold, given comparison of integers. |

458 | { |

459 | any_inverted = true; |

460 | x_tilde[j] = rounded_value; |

461 | |

462 | // the following is simply for statistics. |

463 | if (rounded_value <= new_solution[j]) |

464 | { |

465 | ++num_rounded_down; |

466 | } |

467 | else |

468 | { |

469 | ++num_rounded_up; |

470 | } |

471 | |

472 | // spit out the actual rounding information for debug/analysis purposes. |

473 | if (verbosity(50) == true) |

474 | { |

475 | std::cout << "Rounded value=" << new_solution[j] << " to integral value=" << rounded_value << std::endl; |

476 | } |

477 | } |

478 | } |

479 | } |

480 | } |

481 | |

482 | return any_inverted; |

483 | } |

484 | |

485 | bool feasibilityPump::randomizedRounding(unsigned int & num_rounded_up, |

486 | unsigned int & num_rounded_down, |

487 | int num_vars, |

488 | const double * new_solution, |

489 | double * x_tilde)const |

490 | { |

491 | bool any_inverted = false; |

492 | num_rounded_down = 0; |

493 | num_rounded_up = 0; |

494 | |

495 | RNG * rng = pebbl::gRandomRNG(); |

496 | Uniform sampler(rng, 0.0, 1.0); |

497 | |

498 | for (int j = 0; j < num_vars; ++j) |

499 | { |

500 | if (lp->isContinuous(j) == false) |

501 | { |

502 | if (bGlobal->isInteger(new_solution[j]) == false) |

503 | { |

504 | // if you're <= 0.33 or >= 0.67 in the fractional component, |

505 | // perform deterministic rounding. otherweise, round with |

506 | // equal probability in one or the other direction. |

507 | double fractional_component(new_solution[j] - floor(new_solution[j])); |

508 | double rounded_value; |

509 | if (fractional_component < 0.33) |

510 | { |

511 | rounded_value = floor(new_solution[j]); |

512 | } |

513 | else if (fractional_component > 0.67) |

514 | { |

515 | rounded_value = ceil(new_solution[j]); |

516 | } |

517 | else |

518 | { |

519 | if (sampler() <= 0.5) |

520 | { |

521 | rounded_value = floor(new_solution[j]); |

522 | } |

523 | else |

524 | { |

525 | rounded_value = ceil(new_solution[j]); |

526 | } |

527 | } |

528 | |

529 | if (fabs(rounded_value - x_tilde[j]) > 1e-5) // enough of a threshold, given comparison of integers. |

530 | { |

531 | any_inverted = true; |

532 | x_tilde[j] = rounded_value; |

533 | |

534 | // the following is simply for statistics. |

535 | if (fabs(rounded_value - floor(new_solution[j])) < 1e-5) |

536 | { |

537 | ++num_rounded_down; |

538 | } |

539 | else |

540 | { |

541 | ++num_rounded_up; |

542 | } |

543 | |

544 | // spit out the actual rounding information for debug/analysis purposes. |

545 | if (verbosity(50) == true) |

546 | { |

547 | std::cout << "Rounded value=" << new_solution[j] << " to integral value=" << rounded_value << std::endl; |

548 | } |

549 | } |

550 | } |

551 | } |

552 | } |

553 | |

554 | return any_inverted; |

555 | } |

556 | |

557 | bool feasibilityPump::feasibilityPumpCore(OsiSolverInterface * fp_solver_interface, |

558 | const double * initial_solution, |

559 | double initial_solution_lb, |

560 | int allocated_iterations, |

561 | double allocated_runtime, |

562 | int & consumed_iterations, |

563 | double & consumed_runtime, |

564 | double * incumbent_solution) |

565 | { |

566 | assert(allocated_runtime >= 0.0); |

567 | assert(allocated_iterations >= 0); |

568 | |

569 | double start_time = CPUSeconds(); |

570 | int num_vars = fp_solver_interface->getNumCols(); |

571 | int num_int_vars = problem->numBinaryVars(); |

572 | |

573 | consumed_iterations = 0; |

574 | consumed_runtime = 0.0; |

575 | |

576 | if (FPRoundingType == 0) |

577 | { |

578 | DEBUGPR(10, ucout << "Deterministic rounding enabled for feasibility pump" << std::endl << Flush); |

579 | } |

580 | else // (FPRoundingType == 1) |

581 | { |

582 | DEBUGPR(10, ucout << "Randomized rounding enabled for feasibility pump" << std::endl << Flush); |

583 | } |

584 | |

585 | // initialize the hash weights, which are used to detect cycles in the solution vector x tilde. |

586 | std::vector<int> hash_weights(num_vars); |

587 | for (int i = 0;i < num_vars; ++i) |

588 | { |

589 | RNG * rng = pebbl::gRandomRNG(); |

590 | DUniform<int> sampler(rng, 1, std::numeric_limits<int>::max()); |

591 | hash_weights[i] = sampler(); |

592 | } |

593 | |

594 | const int CYCLE_HISTORY_LENGTH = 100; // could be promoted to a parameter at some point - should be linked to "R". |

595 | std::list<double> x_tilde_hash_history; |

596 | |

597 | // the initial x~ is just the input solution. |

598 | double * x_tilde = new double[num_vars]; |

599 | double hash_value(0.0); |

600 | for (int i = 0; i < num_vars; ++i) |

601 | { |

602 | x_tilde[i] = initial_solution[i]; |

603 | hash_value += x_tilde[i] * hash_weights[i]; |

604 | } |

605 | |

606 | x_tilde_hash_history.push_back(hash_value); |

607 | |

608 | // perform the specified number of feasibility pump iterations. |

609 | bool incumbent_found(false); |

610 | bool time_limited(false); // upon exit, is it due to time limiation? |

611 | bool iteration_limited(false); // upon exit, is it due to iteration limitation? |

612 | int iteration_feasible_found(std::numeric_limits<int>::min()); |

613 | |

614 | // the following is for reporting purposes only, and won't be triggered |

615 | // without some real weird inputs. |

616 | if (allocated_iterations == 0) |

617 | { |

618 | iteration_limited = true; |

619 | } |

620 | |

621 | // in an attempt to accelerate solves... |

622 | CoinWarmStart * warm_start = 0; |

623 | |

624 | double lp_relaxation_objective(fp_solver_interface->getInfinity() * fp_solver_interface->getObjSense()); |

625 | double best_incumbent_objective(fp_solver_interface->getInfinity() * fp_solver_interface->getObjSense()); |

626 | |

627 | for (int current_iteration = 1; current_iteration <= allocated_iterations; ++current_iteration) |

628 | { |

629 | DEBUGPR(100, ucout << "Executing feasibility pump iteration=" << std::setw(5) << current_iteration << std::endl << Flush); |

630 | |

631 | ++consumed_iterations; |

632 | |

633 | DEBUGPR(100, prettyPrint("Reference x~ solution", x_tilde, num_vars, ucout)); |

634 | |

635 | // form the modified objective function for the LP, based on x~ |

636 | unsigned base_offset(0); |

637 | |

638 | for (int j = 0; j < num_vars; ++j) |

639 | { |

640 | // again, the lp instance has the integer information. |

641 | if (lp->isContinuous(j) == true) |

642 | { |

643 | fp_solver_interface->setObjCoeff(j, 0.0); |

644 | } |

645 | else // for now, we assume all integers are binary |

646 | { |

647 | if (x_tilde[j] == 0.0) // the rounded variables should be exactly equal to 0 (which is representable), so tolerances shouldn't be an issue. |

648 | { |

649 | fp_solver_interface->setObjCoeff(j, 1.0); |

650 | } |

651 | else |

652 | { |

653 | base_offset++; |

654 | fp_solver_interface->setObjCoeff(j, -1.0); |

655 | } |

656 | } |

657 | } |

658 | |

659 | const double * coeffs = fp_solver_interface->getObjCoefficients(); |

660 | DEBUGPR(100, prettyPrint("New objective function coefficients in lp*", coeffs, num_vars, ucout)); |

661 | |

662 | // solve LP* |

663 | DEBUGPR(100, ucout << "Solving LP*" << std::endl << Flush); |

664 | |

665 | // NOTE: I had originally tried an "initialSolve()" call on the first iteration, |

666 | // and "resolve()" calls on subsequent iterations. this yielded some really |

667 | // strange behaviors on some instances, including mas74 and opt1217. there |

668 | // may be some confusion over what the term "problem modification" means |

669 | // in OSI land, as the limited documentation indicates you can call "resolve" |

670 | // after a problem modification. |

671 | |

672 | if (warm_start != 0) |

673 | { |

674 | if (fp_solver_interface->setWarmStart(warm_start) == false) |

675 | { |

676 | DEBUGPR(10, ucout << "WARNING: COIN rejected the warm start object" << std::endl << Flush); |

677 | } |

678 | } |

679 | |

680 | fp_solver_interface->messageHandler()->setLogLevel(0); |

681 | fp_solver_interface->initialSolve(); |

682 | |

683 | // save for the next round. |

684 | delete warm_start; |

685 | warm_start = fp_solver_interface->getWarmStart(); // ownership is trasnferred |

686 | |

687 | DEBUGPR(100, ucout << "LP* objective value=" << fp_solver_interface->getObjValue() << std::endl << Flush); |

688 | DEBUGPR(100, ucout << "LP* base offset=" << base_offset << std::endl << Flush); |

689 | const double * new_solution = fp_solver_interface->getColSolution(); |

690 | DEBUGPR(100, prettyPrint("LP* solution", new_solution, num_vars, ucout)); |

691 | |

692 | // compute the L1-norm distance between the new x* and the old x~ |

693 | double dist(0.0); |

694 | for (int j = 0; j < num_vars; ++j) |

695 | { |

696 | if (lp->isContinuous(j) == false) |

697 | { |

698 | dist += fabs(x_tilde[j]-new_solution[j]); |

699 | } |

700 | } |

701 | |

702 | DEBUGPR(10, ucout << "Iteration=" << std::setw(5) << current_iteration << ": Distance=" << std::setw(10) << dist << " " << ": Simplex iterations=" << fp_solver_interface->getIterationCount() << Flush); |

703 | |

704 | // NOTE: The solution can still be integer-feasible, although the |

705 | // distance is > 0. weird, but true. due to tolerances. |

706 | |

707 | // see if the resulting solution is integer. |

708 | bool integer_feasible(true); // until proven otherwise |

709 | for (int j = 0; (integer_feasible == true) && (j < num_vars); ++j) |

710 | { |

711 | if (lp->isContinuous(j) == false) |

712 | { |

713 | // integer tolerances are found in the pebblBase |

714 | // class, as is this static utility method. |

715 | if (bGlobal->isInteger(new_solution[j]) == false) |

716 | { |

717 | integer_feasible = false; |

718 | } |

719 | } |

720 | } |

721 | |

722 | // if the solution is integer feasible, then update the best incumbent |

723 | // and update the objective cut accordingly (assuming continuation is enabled). |

724 | if (integer_feasible == true) |

725 | { |

726 | incumbent_found = true; |

727 | iteration_feasible_found = current_iteration; |

728 | |

729 | // NOTE: this was the original - now keep going after modifying the objective function! |

730 | // break; |

731 | |

732 | const double * new_solution = fp_solver_interface->getColSolution(); |

733 | |

734 | DoubleVector heuristic_solution(num_vars); |

735 | for (int j = 0; j < num_vars; ++j) |

736 | { |

737 | heuristic_solution[j] = new_solution[j]; |

738 | incumbent_solution[j] = new_solution[j]; // save the best-so-far. |

739 | } |

740 | |

741 | double solution_objective = lp->evalVector(heuristic_solution); |

742 | best_incumbent_objective = solution_objective; |

743 | |

744 | DEBUGPR(10, ucout << ", >>>>solution is integer feasible - objective=" << solution_objective << std::endl << Flush); |

745 | |

746 | x_tilde_hash_history.clear(); |

747 | |

748 | // only continue search if the improve-incumbent flag is enabled - |

749 | // otherwise, bail once you have a feasible incumbent. |

750 | if (FPImproveIncumbent == false) |

751 | { |

752 | DEBUGPR(10, ucout << "Incumbent improvement is disabled - terminating search" << std::endl << Flush); |

753 | break; |

754 | } |

755 | |

756 | // TBD - still need to mess with the objective sense - it's currently hardwired for minimization. |

757 | double absolute_gap = solution_objective - initial_solution_lb; |

758 | double new_objective = solution_objective - FPImprovementAlpha * absolute_gap; |

759 | int num_rows = fp_solver_interface->getNumRows(); |

760 | fp_solver_interface->setRowUpper(num_rows-1, new_objective); |

761 | |

762 | DEBUGPR(10, ucout << "Objective upper bound reset to " << new_objective << std::endl << Flush); |

763 | } |

764 | else |

765 | { |

766 | // see if there is a cycle in the x tilde values. |

767 | bool cycle_detected(false); |

768 | unsigned int cycle_length(0); |

769 | std::list<double>::reverse_iterator history_iterator(x_tilde_hash_history.rbegin()); |

770 | double current_hash_value(*history_iterator); |

771 | for (++history_iterator; history_iterator != x_tilde_hash_history.rend(); ++history_iterator) |

772 | { |

773 | const int MINIMUM_CYCLE_LENGTH = 5; // random moves can get you out of simple cases, so let it go a bit. |

774 | ++cycle_length; |

775 | if ((cycle_length >= MINIMUM_CYCLE_LENGTH) && |

776 | (fabs((*history_iterator) - current_hash_value) < 1e-5)) |

777 | { |

778 | DEBUGPR(10, ucout << ", Cycle detected in x tilde vector, cycle length=" << cycle_length << Flush); |

779 | cycle_detected = true; |

780 | break; |

781 | } |

782 | } |

783 | |

784 | // NOTE: there is a hack described in the original feasibility pump |

785 | // paper, in which a random move is applied every FPPerturbationInterval |

786 | // iterations to generate x~, instead of using x*. |

787 | |

788 | if (((current_iteration % FPPerturbationInterval) == 0) || |

789 | (cycle_detected == true)) |

790 | { |

791 | DEBUGPR(10, ucout << ", Applying perturbation" << std::endl << Flush); |

792 | |

793 | RNG * rng = pebbl::gRandomRNG(); |

794 | Uniform sampler(rng, -0.3, 0.7); |

795 | |

796 | for (int j = 0; j < num_vars; ++j) |

797 | { |

798 | if (lp->isContinuous(j) == false) |

799 | { |

800 | double sample = sampler(); |

801 | double result = fabs(new_solution[j] - x_tilde[j]) + std::max(sample,0.0); |

802 | |

803 | if (result >= 0.5) |

804 | { |

805 | x_tilde[j] = (1.0 - x_tilde[j]); |

806 | } |

807 | } |

808 | } |

809 | |

810 | x_tilde_hash_history.clear(); |

811 | } |

812 | else |

813 | { |

814 | // construct the next x tilde via any number of mechanisms, to see if rounding |

815 | // of any x* variables yields new values for the corresponding x~ variables. |

816 | |

817 | bool any_inverted(false); |

818 | unsigned int num_rounded_up(0); |

819 | unsigned int num_rounded_down(0); |

820 | |

821 | if (FPRoundingType == 0) |

822 | { |

823 | any_inverted = deterministicRounding(num_rounded_up, |

824 | num_rounded_down, |

825 | num_vars, |

826 | new_solution, |

827 | x_tilde); |

828 | } |

829 | else // (FPRoundingType == 1) |

830 | { |

831 | any_inverted = randomizedRounding(num_rounded_up, |

832 | num_rounded_down, |

833 | num_vars, |

834 | new_solution, |

835 | x_tilde); |

836 | } |

837 | |

838 | if (any_inverted == true) |

839 | { |

840 | DEBUGPR(10, ucout << ", Variables were inverted - " << num_rounded_down << " rounded down, " << num_rounded_up << " rounded up" << std::endl << Flush); |

841 | } |

842 | else |

843 | { |

844 | // if the new x* rounding failed to change x~, |

845 | // flip some number of variables with the highest |

846 | // value of abs(x* - x~). |

847 | DEBUGPR(100, ucout << "Mean flip count=" << FPMeanFlipCount << std::endl << Flush); |

848 | |

849 | // compute |x* - x~| for each variable and sort in descending order. |

850 | std::vector<std::pair<double,int> > distances; |

851 | distances.reserve(num_vars); |

852 | for (int j = 0; j < num_vars; ++j) |

853 | { |

854 | if (lp->isContinuous(j) == false) |

855 | { |

856 | double this_distance = fabs(new_solution[j]-x_tilde[j]); |

857 | distances.push_back(std::make_pair(this_distance,j)); |

858 | } |

859 | } |

860 | |

861 | std::sort(distances.begin(),distances.end(),compareDists); |

862 | |

863 | // std::cout << "Distance/flip info: "; |

864 | // for (size_t i = 0; i < distances.size(); ++i) |

865 | // { |

866 | // std::cout << "(" << distances[i].first << "," << distances[i].second << ") "; |

867 | // } |

868 | // std::cout << std::endl; |

869 | |

870 | // flip the rand(fpFlipCount/2,3*fpFlipCount/2) variables with the largest distance |

871 | RNG * rng = pebbl::gRandomRNG(); |

872 | DUniform<int> sampler(rng, |

873 | std::max(int(floor(double(FPMeanFlipCount)/2.0)),1), |

874 | int(ceil(3.0*double(FPMeanFlipCount)/2.0))); |

875 | |

876 | int num_to_flip = sampler(); |

877 | // can't flip more variables than there are! |

878 | if (num_to_flip > num_int_vars) |

879 | { |

880 | num_to_flip = num_int_vars; |

881 | } |

882 | |

883 | DEBUGPR(100, ucout << "Actual number of flips to perform=" << num_to_flip << std::endl << Flush); |

884 | |

885 | // normalize flip count to the maximal number of binary variables. |

886 | // NOTE: a lot of these are probably the same, and we should randomize within the blocks. |

887 | num_to_flip = std::min(num_to_flip,int(distances.size())); |

888 | |

889 | for (int j = 0; j < num_to_flip; ++j) |

890 | { |

891 | int this_var = distances[j].second; |

892 | // std::cout << "This distance=" << distances[j].first << std::endl; |

893 | x_tilde[this_var] = 1.0 - x_tilde[this_var]; |

894 | } |

895 | |

896 | DEBUGPR(10, ucout << ", Randomly flipped " << num_to_flip << " variables" << std::endl << Flush); |

897 | |

898 | DEBUGPR(100, prettyPrint("X~ solution after randomized flipping", x_tilde, num_vars, ucout)); |

899 | } |

900 | } |

901 | |

902 | // update x tilde history for cycle detection. |

903 | hash_value = 0.0; |

904 | for (int i = 0; i < num_vars; ++i) |

905 | { |

906 | hash_value += x_tilde[i] * hash_weights[i]; |

907 | } |

908 | |

909 | x_tilde_hash_history.push_back(hash_value); |

910 | |

911 | while (int(x_tilde_hash_history.size()) > CYCLE_HISTORY_LENGTH) |

912 | { |

913 | x_tilde_hash_history.pop_front(); |

914 | } |

915 | } |

916 | |

917 | double elapsed_time = CPUSeconds() - start_time; |

918 | if (elapsed_time >= allocated_runtime) |

919 | { |

920 | time_limited = true; |

921 | |

922 | break; |

923 | } |

924 | |

925 | if (current_iteration == allocated_iterations) |

926 | { |

927 | DEBUGPR(10, ucout << "Feasibility pump terminating after exceeeding allocated number of iterations" << std::endl << Flush); |

928 | |

929 | iteration_limited = true; |

930 | // the loop will terminate naturally at this point - no break is necessary. |

931 | } |

932 | |

933 | DEBUGPR(50, ucout << std::endl << Flush); |

934 | } |

935 | |

936 | double end_time = CPUSeconds(); |

937 | |

938 | if (incumbent_found == true) |

939 | { |

940 | DoubleVector solution_to_evaluate(num_vars); |

941 | for (int j = 0; j < num_vars; ++j) |

942 | { |

943 | solution_to_evaluate[j] = incumbent_solution[j]; |

944 | } |

945 | |

946 | double solution_objective = lp->evalVector(solution_to_evaluate); |

947 | |

948 | DEBUGPR(10, ucout << "Feasibility pump found feasible integer solution: Iterations=" << iteration_feasible_found << ", "); |

949 | DEBUGPR(10, ucout << "Objective value=" << solution_objective << ", "); |

950 | DEBUGPR(10, ucout << "CPU time=" << end_time - start_time << " seconds" << std::endl << Flush); |

951 | } |

952 | else |

953 | { |

954 | DEBUGPR(10, ucout << std::endl << "Feasibility pump failed to find feasible integer solution: "); |

955 | if (iteration_limited == true) |

956 | { |

957 | DEBUGPR(10, ucout << "Iteration limit=" << allocated_iterations << " exceeded"); |

958 | } |

959 | else // time_limited == true |

960 | { |

961 | DEBUGPR(10, ucout << "Time limit=" << allocated_runtime << " seconds exceeded"); |

962 | } |

963 | DEBUGPR(10, ucout << std::endl << Flush); |

964 | } |

965 | |

966 | delete warm_start; |

967 | |

968 | delete [] x_tilde; |

969 | |

970 | consumed_runtime = CPUSeconds() - start_time; |

971 | |

972 | return (incumbent_found == true); |

973 | } |

974 | |

975 | |

976 | } // namespace pico |

**Note:**See TracBrowser for help on using the repository browser.