1

2

3

4

5

6

7

8

9

10

11

12

13

14

15

16

17

18

19

20

21

22

23

24

25

26

27

28

29

30

31

32

33

34

35

36

37

38

39

40

41

42

43

44

45

46

47

48

49

50

51

52

53

54

55

56

57

58

59

60

61

62

63

64

65

66

67

68

69

70

71

72

73

74

75

76

77

78

79

80

81

82

83

84

85

86

87

88

89

90

91

92

93

94

95

96

97

98

99

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121#ifndef ORTOOLS_GRAPH_ONE_TREE_LOWER_BOUND_H_

122#define ORTOOLS_GRAPH_ONE_TREE_LOWER_BOUND_H_

123

124#include <cmath>

125#include <cstdint>

126#include <limits>

127#include <set>

128#include <utility>

129#include <vector>

130

131#include "absl/log/log.h"

132#include "absl/types/span.h"

136

138

139

140

141

142

143

144

145

146

147

148

149

150

151

152

153

154

155

156template <typename CostType>

157class VolgenantJonkerEvaluator {

158 public:

159 VolgenantJonkerEvaluator(int number_of_nodes, int max_iterations)

160 : step1_initialized_(false),

161 step1_(0),

162 iteration_(0),

164 : MaxIterations(number_of_nodes)),

165 number_of_nodes_(number_of_nodes) {}

166

167 bool Next() { return iteration_++ < max_iterations_; }

168

178

186

187 void OnNewWMax(CostType one_tree_cost) { UpdateStep(one_tree_cost); }

188

189 private:

190

191

192 static int MaxIterations(int number_of_nodes) {

193 return static_cast<int>(28 * std::pow(number_of_nodes, 0.62));

194 }

195

196 void UpdateStep(CostType one_tree_cost) {

197 step1_ = one_tree_cost / (2 * number_of_nodes_);

198 }

199

200 bool step1_initialized_;

201 double step1_;

202 int iteration_;

203 const int max_iterations_;

204 const int number_of_nodes_;

205};

206

207

208

209template <typename CostType, typename CostFunction>

210class HeldWolfeCrowderEvaluator {

211 public:

212 HeldWolfeCrowderEvaluator(int number_of_nodes, const CostFunction& cost)

213 : iteration_(0),

214 number_of_iterations_(2 * number_of_nodes),

215 upper_bound_(0),

216 lambda_(2.0),

217 step_(0) {

218

219

220 ChristofidesPathSolver<CostType, int64_t, int, CostFunction> solver(

221 number_of_nodes, cost);

222 upper_bound_ = solver.TravelingSalesmanCost().value();

223 }

224

225 bool Next() {

226 const int min_iterations = 2;

227 if (iteration_ >= number_of_iterations_) {

228 number_of_iterations_ /= 2;

229 if (number_of_iterations_ < min_iterations) return false;

230 iteration_ = 0;

231 lambda_ /= 2;

232 } else {

233 ++iteration_;

234 }

235 return true;

236 }

237

238 double GetStep() const { return step_; }

239

240 void OnOneTree(CostType one_tree_cost, double w,

241 absl::Span<const int> degrees) {

242 double norm = 0;

243 for (int degree : degrees) {

244 const double delta = degree - 2;

245 norm += delta * delta;

246 }

247 step_ = lambda_ * (upper_bound_ - w) / norm;

248 }

249

250 void OnNewWMax(CostType one_tree_cost) {}

251

252 private:

253 int iteration_;

254 int number_of_iterations_;

255 CostType upper_bound_;

256 double lambda_;

257 double step_;

258};

259

260

261

262

263

264

265template <typename CostFunction>

266std::set<std::pair<int, int>> NearestNeighbors(int number_of_nodes,

267 int number_of_neighbors,

268 const CostFunction& cost) {

269 using CostType = decltype(cost(0, 0));

270 std::set<std::pair<int, int>> nearest;

271 for (int i = 0; i < number_of_nodes; ++i) {

272 std::vector<std::pair<CostType, int>> neighbors;

273 neighbors.reserve(number_of_nodes - 1);

274 for (int j = 0; j < number_of_nodes; ++j) {

275 if (i != j) {

276 neighbors.emplace_back(cost(i, j), j);

277 }

278 }

279 int size = neighbors.size();

280 if (number_of_neighbors < size) {

281 std::nth_element(neighbors.begin(),

282 neighbors.begin() + number_of_neighbors - 1,

283 neighbors.end());

284 size = number_of_neighbors;

285 }

286 for (int j = 0; j < size; ++j) {

287 nearest.insert({i, neighbors[j].second});

288 nearest.insert({neighbors[j].second, i});

289 }

290 }

291 return nearest;

292}

293

294

295

296template <typename CostFunction>

297void AddArcsFromMinimumSpanningTree(int number_of_nodes,

298 const CostFunction& cost,

299 std::set<std::pair<int, int>>* arcs) {

300 util::CompleteGraph<int, int> graph(number_of_nodes);

301 const std::vector<int> mst =

303 return cost(graph.Tail(arc), graph.Head(arc));

304 });

305 for (int arc : mst) {

306 arcs->insert({graph.Tail(arc), graph.Head(arc)});

307 arcs->insert({graph.Head(arc), graph.Tail(arc)});

308 }

309}

310

311

312

313template <typename CostFunction, typename GraphType, typename AcceptFunction>

314int GetNodeMinimizingEdgeCostToSource(const GraphType& graph, int source,

315 const CostFunction& cost,

316 AcceptFunction accept) {

317 int best_node = -1;

318 double best_edge_cost = 0;

319 for (const auto node : graph.AllNodes()) {

320 if (accept(node)) {

321 const double edge_cost = cost(node, source);

322 if (best_node == -1 || edge_cost < best_edge_cost) {

323 best_node = node;

324 best_edge_cost = edge_cost;

325 }

326 }

327 }

328 return best_node;

329}

330

331

332

333

334template <typename CostFunction, typename GraphType, typename CostType>

335std::vector<int> ComputeOneTree(const GraphType& graph,

336 const CostFunction& cost,

337 absl::Span<const double> weights,

338 absl::Span<const int> sorted_arcs,

339 CostType* one_tree_cost) {

340 const auto weighed_cost = [&cost, weights](int from, int to) {

341 return cost(from, to) + weights[from] + weights[to];

342 };

343

344 std::vector<int> mst;

345 if (!sorted_arcs.empty()) {

347 sorted_arcs);

348 } else {

350 graph, [&weighed_cost, &graph](int arc) {

351 return weighed_cost(graph.Tail(arc), graph.Head(arc));

352 });

353 }

354 std::vector<int> degrees(graph.num_nodes() + 1, 0);

355 *one_tree_cost = 0;

356 for (int arc : mst) {

357 degrees[graph.Head(arc)]++;

358 degrees[graph.Tail(arc)]++;

359 *one_tree_cost += cost(graph.Tail(arc), graph.Head(arc));

360 }

361

362

363 const int extra_node = graph.num_nodes();

364 const auto update_one_tree = [extra_node, one_tree_cost, &degrees,

365 &cost](int node) {

366 *one_tree_cost += cost(node, extra_node);

369 };

370 const int node = GetNodeMinimizingEdgeCostToSource(

371 graph, extra_node, weighed_cost,

372 [extra_node](int n) { return n != extra_node; });

373 update_one_tree(node);

374 update_one_tree(GetNodeMinimizingEdgeCostToSource(

375 graph, extra_node, weighed_cost,

376 [extra_node, node](int n) { return n != extra_node && n != node; }));

378}

379

380

381template <typename CostFunction, typename Algorithm>

382double ComputeOneTreeLowerBoundWithAlgorithm(int number_of_nodes,

383 int nearest_neighbors,

384 const CostFunction& cost,

385 Algorithm* algorithm) {

386 if (number_of_nodes < 2) return 0;

387 if (number_of_nodes == 2) return cost(0, 1) + cost(1, 0);

388 using CostType = decltype(cost(0, 0));

389 auto nearest = NearestNeighbors(number_of_nodes - 1, nearest_neighbors, cost);

390

391

392

393 AddArcsFromMinimumSpanningTree(number_of_nodes - 1, cost, &nearest);

394 util::ListGraph<int, int> graph(number_of_nodes - 1, nearest.size());

395 for (const auto& arc : nearest) {

396 graph.AddArc(arc.first, arc.second);

397 }

398 std::vector<double> weights(number_of_nodes, 0);

399 std::vector<double> best_weights(number_of_nodes, 0);

400 double max_w = -std::numeric_limits<double>::infinity();

401 double w = 0;

402

403 while (algorithm->Next()) {

404 CostType one_tree_cost = 0;

405 const std::vector<int> degrees =

406 ComputeOneTree(graph, cost, weights, {}, &one_tree_cost);

407 algorithm->OnOneTree(one_tree_cost, w, degrees);

408 w = one_tree_cost;

409 for (int j = 0; j < number_of_nodes; ++j) {

410 w += weights[j] * (degrees[j] - 2);

411 }

412 if (w > max_w) {

413 max_w = w;

414 best_weights = weights;

415 algorithm->OnNewWMax(one_tree_cost);

416 }

417 const double step = algorithm->GetStep();

418 for (int j = 0; j < number_of_nodes; ++j) {

419 weights[j] += step * (degrees[j] - 2);

420 }

421 }

422

423

424

425 util::CompleteGraph<int, int> complete_graph(number_of_nodes - 1);

426 CostType one_tree_cost = 0;

427

428

429

430 const std::vector<int> degrees =

431 ComputeOneTree(complete_graph, cost, best_weights, {}, &one_tree_cost);

432 w = one_tree_cost;

433 for (int j = 0; j < number_of_nodes; ++j) {

434 w += best_weights[j] * (degrees[j] - 2);

435 }

436 return w;

437}

438

439

440struct TravelingSalesmanLowerBoundParameters {

441 enum Algorithm {

442 VolgenantJonker,

443 HeldWolfeCrowder,

444 };

445

446 Algorithm algorithm = VolgenantJonker;

447

448

449 int volgenant_jonker_iterations = 0;

450

451 int nearest_neighbors = 40;

452};

453

454

455template <typename CostFunction>

456double ComputeOneTreeLowerBoundWithParameters(

457 int number_of_nodes, const CostFunction& cost,

458 const TravelingSalesmanLowerBoundParameters& parameters) {

459 using CostType = decltype(cost(0, 0));

460 switch (parameters.algorithm) {

461 case TravelingSalesmanLowerBoundParameters::VolgenantJonker: {

462 VolgenantJonkerEvaluator<CostType> algorithm(

463 number_of_nodes, parameters.volgenant_jonker_iterations);

464 return ComputeOneTreeLowerBoundWithAlgorithm(

465 number_of_nodes, parameters.nearest_neighbors, cost, &algorithm);

466 break;

467 }

468 case TravelingSalesmanLowerBoundParameters::HeldWolfeCrowder: {

469 HeldWolfeCrowderEvaluator<CostType, CostFunction> algorithm(

470 number_of_nodes, cost);

471 return ComputeOneTreeLowerBoundWithAlgorithm(

472 number_of_nodes, parameters.nearest_neighbors, cost, &algorithm);

473 }

474 default:

475 LOG(ERROR) << "Unsupported algorithm: " << parameters.algorithm;

476 return 0;

477 }

478}

479

480

481

482

483template <typename CostFunction>

484double ComputeOneTreeLowerBound(int number_of_nodes, const CostFunction& cost) {

485 TravelingSalesmanLowerBoundParameters parameters;

486 return ComputeOneTreeLowerBoundWithParameters(number_of_nodes, cost,

487 parameters);

488}

489

490}

491

492#endif