57 public:

66

67

68

69

70

71

72

76

77

78

80

81

82

84

85

86 absl::Status Solve();

87

88 private:

89

90 template <typename T>

91 T SafeAdd(T a, T b) {

92

93 if constexpr (std::is_same_v<std::remove_cv_t<std::remove_reference_t<T>>,

94 int64_t> == true) {

96 } else {

97 return a + b;

98 }

99 }

100

101

103

104

105 CompleteGraph<NodeIndex, ArcIndex> graph_;

106

107

108 const CostFunction costs_;

109

110

111 CostType tsp_cost_;

112

113

114 std::vector<NodeIndex> tsp_path_;

115

116

117 bool solved_;

118};

125 const WeightFunctionType& weight) {

126 using ArcIndex = typename GraphType::ArcIndex;

127 using NodeIndex = typename GraphType::NodeIndex;

128 if constexpr (!std::is_integral_v<CostType>) {

129 DLOG(WARNING) << "Weights are being cast to int64_t. This might result "

130 "in loss of precision or overflows.";

131 }

133 for (NodeIndex tail : graph.AllNodes()) {

134 for (const ArcIndex arc : graph.OutgoingArcs(tail)) {

135 const NodeIndex head = graph.Head(arc);

136

137 if (tail < head) {

139 }

140 }

141 }

145 return absl::InvalidArgumentError(

146 absl::StrCat("Perfect matching failed: ", status));

147 }

148 std::vector<std::pair<NodeIndex, NodeIndex>> match;

149 for (NodeIndex tail : graph.AllNodes()) {

151 if (tail < head) {

152 match.emplace_back(tail, head);

153 }

154 }

155 return match;

156}

167 const WeightFunctionType& weight) {

168 using ArcIndex = typename GraphType::ArcIndex;

169 using NodeIndex = typename GraphType::NodeIndex;

172

173

174

175

176 std::vector<int> variable_indices(graph.num_arcs(), -1);

177 for (NodeIndex node : graph.AllNodes()) {

178

179 for (const ArcIndex arc : graph.OutgoingArcs(node)) {

180 const NodeIndex head = graph.Head(arc);

181 if (node < head) {

188 }

189 }

190

191

195 }

196 for (NodeIndex node : graph.AllNodes()) {

197 for (const ArcIndex arc : graph.OutgoingArcs(node)) {

198 const NodeIndex head = graph.Head(arc);

199 if (node < head) {

200 const int arc_var = variable_indices[arc];

201 DCHECK_GE(arc_var, 0);

208 }

209 }

210 }

211#if defined(USE_SCIP)

212 MPSolver mp_solver("MatchingWithSCIP",

214#elif defined(USE_CBC)

215 MPSolver mp_solver("MatchingWithCBC",

217#endif

218 std::string error;

222 return absl::InvalidArgumentError("MIP-based matching failed");

223 }

226 std::vector<std::pair<NodeIndex, NodeIndex>> matching;

227 for (ArcIndex arc = 0; arc < variable_indices.size(); ++arc) {

228 const int arc_var = variable_indices[arc];

229 if (arc_var >= 0 && response.variable_value(arc_var) > .9) {

230 DCHECK_GE(response.variable_value(arc_var), 1.0 - 1e-4);

231 matching.emplace_back(graph.Tail(arc), graph.Head(arc));

232 }

233 }

234 return matching;

235}

243 graph_(num_nodes),

244 costs_(std::move(costs)),

245 tsp_cost_(0),

246 solved_(false) {}

252 if (!solved_) {

253 const absl::Status status = Solve();

254 if (!status.ok()) return status;

255 }

256 return tsp_cost_;

257}

263 if (!solved_) {

264 const absl::Status status = Solve();

265 if (!status.ok()) return status;

266 }

267 return tsp_path_;

268}

274 const NodeIndex num_nodes = graph_.num_nodes();

275 tsp_path_.clear();

276 tsp_cost_ = 0;

277 if (num_nodes == 1) {

278 tsp_path_ = {0, 0};

279 }

280 if (num_nodes <= 1) {

281 return absl::OkStatus();

282 }

283

284 const std::vector<ArcIndex> mst =

286 return costs_(graph_.Tail(arc), graph_.Head(arc));

287 });

288

289 std::vector<NodeIndex> degrees(num_nodes, 0);

290 for (ArcIndex arc : mst) {

291 degrees[graph_.Tail(arc)]++;

292 degrees[graph_.Head(arc)]++;

293 }

294 std::vector<NodeIndex> odd_degree_nodes;

295 for (int i = 0; i < degrees.size(); ++i) {

296 if (degrees[i] % 2 != 0) {

297 odd_degree_nodes.push_back(i);

298 }

299 }

300

301 const NodeIndex reduced_size = odd_degree_nodes.size();

302 DCHECK_NE(0, reduced_size);

304 std::vector<std::pair<NodeIndex, NodeIndex>> closure_arcs;

305 switch (matching_) {

308 reduced_graph, [this, &reduced_graph,

310 return costs_(odd_degree_nodes[reduced_graph.Tail(arc)],

311 odd_degree_nodes[reduced_graph.Head(arc)]);

312 });

313 if (!result.ok()) {

314 return result.status();

315 }

316 result->swap(closure_arcs);

317 break;

318 }

319#if defined(USE_CBC) || defined(USE_SCIP)

322 reduced_graph, [this, &reduced_graph,

324 return costs_(odd_degree_nodes[reduced_graph.Tail(arc)],

325 odd_degree_nodes[reduced_graph.Head(arc)]);

326 });

327 if (!result.ok()) {

328 return result.status();

329 }

330 result->swap(closure_arcs);

331 break;

332 }

333#endif

335

336

337 std::vector<ArcIndex> ordered_arcs(reduced_graph.num_arcs());

338 std::vector<CostType> ordered_arc_costs(reduced_graph.num_arcs(), 0);

339 for (const ArcIndex arc : reduced_graph.AllForwardArcs()) {

340 ordered_arcs[arc] = arc;

341 ordered_arc_costs[arc] =

342 costs_(odd_degree_nodes[reduced_graph.Tail(arc)],

343 odd_degree_nodes[reduced_graph.Head(arc)]);

344 }

345 absl::c_sort(ordered_arcs,

346 [&ordered_arc_costs](ArcIndex arc_a, ArcIndex arc_b) {

347 return ordered_arc_costs[arc_a] < ordered_arc_costs[arc_b];

348 });

349 std::vector<bool> touched_nodes(reduced_size, false);

350 for (ArcIndex arc_index = 0; closure_arcs.size() * 2 < reduced_size;

351 ++arc_index) {

352 const ArcIndex arc = ordered_arcs[arc_index];

355 if (head != tail && !touched_nodes[tail] && !touched_nodes[head]) {

356 touched_nodes[tail] = true;

357 touched_nodes[head] = true;

358 closure_arcs.emplace_back(tail, head);

359 }

360 }

361 break;

362 }

363 }

364

365

366

368 num_nodes, closure_arcs.size() + mst.size());

369 for (ArcIndex arc : mst) {

370 egraph.AddArc(graph_.Tail(arc), graph_.Head(arc));

371 }

372 for (const auto arc : closure_arcs) {

373 egraph.AddArc(odd_degree_nodes[arc.first], odd_degree_nodes[arc.second]);

374 }

375 std::vector<bool> touched(num_nodes, false);

378 if (touched[node]) continue;

379 touched[node] = true;

380 tsp_cost_ = SafeAdd(tsp_cost_,

381 tsp_path_.empty() ? 0 : costs_(tsp_path_.back(), node));

382 tsp_path_.push_back(node);

383 }

384 tsp_cost_ =

385 SafeAdd(tsp_cost_, tsp_path_.empty() ? 0 : costs_(tsp_path_.back(), 0));

386 tsp_path_.push_back(0);

387 solved_ = true;

388 return absl::OkStatus();

389}

Definition christofides.h:56

absl::Status Solve()

Definition christofides.h:273

ChristofidesPathSolver(NodeIndex num_nodes, CostFunction costs)

Definition christofides.h:241

void SetMatchingAlgorithm(MatchingAlgorithm matching)

Definition christofides.h:73

absl::StatusOr< CostType > TravelingSalesmanCost()

Definition christofides.h:251

absl::StatusOr< const std::vector< NodeIndex > & > TravelingSalesmanPath()

Definition christofides.h:262

MatchingAlgorithm

Definition christofides.h:58

@ MINIMUM_WEIGHT_MATCHING_WITH_MIP

Definition christofides.h:61

@ MINIMAL_WEIGHT_MATCHING

Definition christofides.h:63

@ MINIMUM_WEIGHT_MATCHING

Definition christofides.h:59

void set_lower_bound(double value)

void add_coefficient(double value)

void set_upper_bound(double value)

void add_var_index(::int32_t value)

int variable_size() const

void set_maximize(bool value)

::operations_research::MPConstraintProto *PROTOBUF_NONNULL add_constraint()

::operations_research::MPConstraintProto *PROTOBUF_NONNULL mutable_constraint(int index)

::operations_research::MPVariableProto *PROTOBUF_NONNULL add_variable()

double variable_value(int index) const

ResultStatus Solve()

Solves the problem using the default parameter values.

MPSolverResponseStatus LoadModelFromProto(const MPModelProto &input_model, std::string *error_message, bool clear_names=true)

@ SCIP_MIXED_INTEGER_PROGRAMMING

@ CBC_MIXED_INTEGER_PROGRAMMING

void FillSolutionResponseProto(MPSolutionResponse *response) const

Encodes the current solution in a solution response protocol buffer.

void set_upper_bound(double value)

void set_is_integer(bool value)

void set_objective_coefficient(double value)

void set_lower_bound(double value)

int Match(int node) const

void AddEdgeWithCost(int tail, int head, int64_t cost)

ABSL_MUST_USE_RESULT Status Solve()

IntegerRange< ArcIndex > AllForwardArcs() const

ArcIndexType num_arcs() const

NodeIndexType Tail(ArcIndexType arc) const

NodeIndexType Head(ArcIndexType arc) const

ArcIndexType AddArc(NodeIndexType tail, NodeIndexType head)

bool IsEulerianGraph(const Graph &graph, bool assume_connectivity=true)

int64_t CapAdd(int64_t x, int64_t y)

absl::StatusOr< std::vector< std::pair< typename GraphType::NodeIndex, typename GraphType::NodeIndex > > > ComputeMinimumWeightMatchingWithMIP(const GraphType &graph, const WeightFunctionType &weight)

Definition christofides.h:166

std::vector< NodeIndex > BuildEulerianTourFromNode(const Graph &graph, NodeIndex root, bool assume_connectivity=true)

BlossomGraph::NodeIndex NodeIndex

std::vector< typename Graph::ArcIndex > BuildPrimMinimumSpanningTree(const Graph &graph, const ArcValue &arc_value)

absl::StatusOr< std::vector< std::pair< typename GraphType::NodeIndex, typename GraphType::NodeIndex > > > ComputeMinimumWeightMatching(const GraphType &graph, const WeightFunctionType &weight)

Definition christofides.h:124

trees with all degrees equal w the current value of degrees