C++ Reference

C++ Reference: Graph

christofides.h
Go to the documentation of this file.
1// Copyright 2010-2018 Google LLC
2// Licensed under the Apache License, Version 2.0 (the "License");
3// you may not use this file except in compliance with the License.
4// You may obtain a copy of the License at
5//
6// http://www.apache.org/licenses/LICENSE-2.0
7//
8// Unless required by applicable law or agreed to in writing, software
9// distributed under the License is distributed on an "AS IS" BASIS,
10// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
11// See the License for the specific language governing permissions and
12// limitations under the License.
13
14// ChristofidesPathSolver computes an approximate solution to the Traveling
15// Salesman Problen using the Christofides algorithm (c.f.
16// https://en.wikipedia.org/wiki/Christofides_algorithm).
17// Note that the algorithm guarantees finding a solution within 3/2 of the
18// optimum. Its complexity is O(n^2 * log(n)) where n is the number of nodes.
19
20#ifndef OR_TOOLS_GRAPH_CHRISTOFIDES_H_
21#define OR_TOOLS_GRAPH_CHRISTOFIDES_H_
22
23#include "absl/status/status.h"
24#include "absl/status/statusor.h"
25#include "ortools/base/integral_types.h"
26#include "ortools/base/logging.h"
28#include "ortools/graph/graph.h"
30#include "ortools/graph/perfect_matching.h"
31#include "ortools/linear_solver/linear_solver.h"
32#include "ortools/linear_solver/linear_solver.pb.h"
33#include "ortools/util/saturated_arithmetic.h"
34
36
37using ::util::CompleteGraph;
38
39template <typename CostType, typename ArcIndex = int64,
40 typename NodeIndex = int32,
41 typename CostFunction = std::function<CostType(NodeIndex, NodeIndex)>>
43 public:
44 enum class MatchingAlgorithm {
45 MINIMUM_WEIGHT_MATCHING,
46#if defined(USE_CBC) || defined(USE_SCIP)
47 MINIMUM_WEIGHT_MATCHING_WITH_MIP,
48#endif // defined(USE_CBC) || defined(USE_SCIP)
49 MINIMAL_WEIGHT_MATCHING,
50 };
51 ChristofidesPathSolver(NodeIndex num_nodes, CostFunction costs);
52
53 // Sets the matching algorith to use. A minimum weight perfect matching
54 // (MINIMUM_WEIGHT_MATCHING) guarantees the 3/2 upper bound to the optimal
55 // solution. A minimal weight perfect matching (MINIMAL_WEIGHT_MATCHING)
56 // finds a locally minimal weight matching which does not offer any bound
57 // guarantee but, as of 1/2017, is orders of magnitude faster than the
58 // minimum matching.
59 // By default, MINIMAL_WEIGHT_MATCHING is selected.
60 // TODO(user): Change the default when minimum matching gets faster.
62 matching_ = matching;
63 }
64
65 // Returns the cost of the approximate TSP tour.
66 CostType TravelingSalesmanCost();
67
68 // Returns the approximate TSP tour.
69 std::vector<NodeIndex> TravelingSalesmanPath();
70
71 // Runs the Christofides algorithm. Returns true if a solution was found,
72 // false otherwise.
73 bool Solve();
74
75 private:
76 // Safe addition operator to avoid overflows when possible.
77 //template <typename T>
78 //T SafeAdd(T a, T b) {
79 // return a + b;
80 //}
81 //template <>
82 int64 SafeAdd(int64 a, int64 b) {
83 return CapAdd(a, b);
84 }
85
86 // Matching algorithm to use.
87 MatchingAlgorithm matching_;
88
89 // The complete graph on the nodes of the problem.
90 CompleteGraph<NodeIndex, ArcIndex> graph_;
91
92 // Function returning the cost between nodes of the problem.
93 const CostFunction costs_;
94
95 // The cost of the computed TSP path.
96 CostType tsp_cost_;
97
98 // The path of the computed TSP,
99 std::vector<NodeIndex> tsp_path_;
100
101 // True if the TSP has been solved, false otherwise.
102 bool solved_;
103};
104
105// Computes a minimum weight perfect matching on an undirected graph.
106template <typename WeightFunctionType, typename GraphType>
107absl::StatusOr<std::vector<
108 std::pair<typename GraphType::NodeIndex, typename GraphType::NodeIndex>>>
109ComputeMinimumWeightMatching(const GraphType& graph,
110 const WeightFunctionType& weight) {
111 using ArcIndex = typename GraphType::ArcIndex;
112 using NodeIndex = typename GraphType::NodeIndex;
113 MinCostPerfectMatching matching(graph.num_nodes());
114 for (NodeIndex tail : graph.AllNodes()) {
115 for (const ArcIndex arc : graph.OutgoingArcs(tail)) {
116 const NodeIndex head = graph.Head(arc);
117 // Adding both arcs is redudant for MinCostPerfectMatching.
118 if (tail < head) {
119 matching.AddEdgeWithCost(tail, head, weight(arc));
120 }
121 }
122 }
123 MinCostPerfectMatching::Status status = matching.Solve();
124 if (status != MinCostPerfectMatching::OPTIMAL) {
125 return absl::InvalidArgumentError("Perfect matching failed");
126 }
127 std::vector<std::pair<NodeIndex, NodeIndex>> match;
128 for (NodeIndex tail : graph.AllNodes()) {
129 const NodeIndex head = matching.Match(tail);
130 if (tail < head) { // Both arcs are matched for a given edge, we keep one.
131 match.emplace_back(tail, head);
132 }
133 }
134 return match;
135}
136
137#if defined(USE_CBC) || defined(USE_SCIP)
138// Computes a minimum weight perfect matching on an undirected graph using a
139// Mixed Integer Programming model.
140// TODO(user): Handle infeasible cases if this algorithm is used outside of
141// Christofides.
142template <typename WeightFunctionType, typename GraphType>
143absl::StatusOr<std::vector<
144 std::pair<typename GraphType::NodeIndex, typename GraphType::NodeIndex>>>
146 const WeightFunctionType& weight) {
147 using ArcIndex = typename GraphType::ArcIndex;
148 using NodeIndex = typename GraphType::NodeIndex;
149 MPModelProto model;
150 model.set_maximize(false);
151 // The model is composed of Boolean decision variables to select matching arcs
152 // and constraints ensuring that each node appears in exactly one selected
153 // arc. The objective is to minimize the sum of the weights of selected arcs.
154 // It is assumed the graph is symmetrical.
155 std::vector<int> variable_indices(graph.num_arcs(), -1);
156 for (NodeIndex node : graph.AllNodes()) {
157 // Creating arc-selection Boolean variable.
158 for (const ArcIndex arc : graph.OutgoingArcs(node)) {
159 const NodeIndex head = graph.Head(arc);
160 if (node < head) {
161 variable_indices[arc] = model.variable_size();
162 MPVariableProto* const arc_var = model.add_variable();
163 arc_var->set_lower_bound(0);
164 arc_var->set_upper_bound(1);
165 arc_var->set_is_integer(true);
166 arc_var->set_objective_coefficient(weight(arc));
167 }
168 }
169 // Creating matching constraint:
170 // for all node i, sum(j) arc(i,j) + sum(j) arc(j,i) = 1
171 MPConstraintProto* const one_of_ct = model.add_constraint();
172 one_of_ct->set_lower_bound(1);
173 one_of_ct->set_upper_bound(1);
174 }
175 for (NodeIndex node : graph.AllNodes()) {
176 for (const ArcIndex arc : graph.OutgoingArcs(node)) {
177 const NodeIndex head = graph.Head(arc);
178 if (node < head) {
179 const int arc_var = variable_indices[arc];
180 DCHECK_GE(arc_var, 0);
181 MPConstraintProto* one_of_ct = model.mutable_constraint(node);
182 one_of_ct->add_var_index(arc_var);
183 one_of_ct->add_coefficient(1);
184 one_of_ct = model.mutable_constraint(head);
185 one_of_ct->add_var_index(arc_var);
186 one_of_ct->add_coefficient(1);
187 }
188 }
189 }
190#if defined(USE_SCIP)
191 MPSolver mp_solver("MatchingWithSCIP",
192 MPSolver::SCIP_MIXED_INTEGER_PROGRAMMING);
193#elif defined(USE_CBC)
194 MPSolver mp_solver("MatchingWithCBC",
195 MPSolver::CBC_MIXED_INTEGER_PROGRAMMING);
196#endif
197 std::string error;
198 mp_solver.LoadModelFromProto(model, &error);
199 MPSolver::ResultStatus status = mp_solver.Solve();
200 if (status != MPSolver::OPTIMAL) {
201 return absl::InvalidArgumentError("MIP-based matching failed");
202 }
203 MPSolutionResponse response;
204 mp_solver.FillSolutionResponseProto(&response);
205 std::vector<std::pair<NodeIndex, NodeIndex>> matching;
206 for (ArcIndex arc = 0; arc < variable_indices.size(); ++arc) {
207 const int arc_var = variable_indices[arc];
208 if (arc_var >= 0 && response.variable_value(arc_var) > .9) {
209 DCHECK_GE(response.variable_value(arc_var), 1.0 - 1e-4);
210 matching.emplace_back(graph.Tail(arc), graph.Head(arc));
211 }
212 }
213 return matching;
214}
215#endif // defined(USE_CBC) || defined(USE_SCIP)
216
217template <typename CostType, typename ArcIndex, typename NodeIndex,
218 typename CostFunction>
220 ChristofidesPathSolver(NodeIndex num_nodes, CostFunction costs)
221 : matching_(MatchingAlgorithm::MINIMAL_WEIGHT_MATCHING),
222 graph_(num_nodes),
223 costs_(std::move(costs)),
224 tsp_cost_(0),
225 solved_(false) {}
226
227template <typename CostType, typename ArcIndex, typename NodeIndex,
228 typename CostFunction>
229CostType ChristofidesPathSolver<CostType, ArcIndex, NodeIndex,
230 CostFunction>::TravelingSalesmanCost() {
231 if (!solved_) {
232 bool const ok = Solve();
233 DCHECK(ok);
234 }
235 return tsp_cost_;
236}
237
238template <typename CostType, typename ArcIndex, typename NodeIndex,
239 typename CostFunction>
240std::vector<NodeIndex> ChristofidesPathSolver<
241 CostType, ArcIndex, NodeIndex, CostFunction>::TravelingSalesmanPath() {
242 if (!solved_) {
243 const bool ok = Solve();
244 DCHECK(ok);
245 }
246 return tsp_path_;
247}
248
249template <typename CostType, typename ArcIndex, typename NodeIndex,
250 typename CostFunction>
252 CostFunction>::Solve() {
253 const NodeIndex num_nodes = graph_.num_nodes();
254 tsp_path_.clear();
255 tsp_cost_ = 0;
256 if (num_nodes == 1) {
257 tsp_path_ = {0, 0};
258 }
259 if (num_nodes <= 1) {
260 return true;
261 }
262 // Compute Minimum Spanning Tree.
263 const std::vector<ArcIndex> mst =
264 BuildPrimMinimumSpanningTree(graph_, [this](ArcIndex arc) {
265 return costs_(graph_.Tail(arc), graph_.Head(arc));
266 });
267 // Detect odd degree nodes.
268 std::vector<NodeIndex> degrees(num_nodes, 0);
269 for (ArcIndex arc : mst) {
270 degrees[graph_.Tail(arc)]++;
271 degrees[graph_.Head(arc)]++;
272 }
273 std::vector<NodeIndex> odd_degree_nodes;
274 for (int i = 0; i < degrees.size(); ++i) {
275 if (degrees[i] % 2 != 0) {
276 odd_degree_nodes.push_back(i);
277 }
278 }
279 // Find minimum-weight perfect matching on odd-degree-node complete graph.
280 // TODO(user): Make this code available as an independent algorithm.
281 const NodeIndex reduced_size = odd_degree_nodes.size();
282 DCHECK_NE(0, reduced_size);
283 CompleteGraph<NodeIndex, ArcIndex> reduced_graph(reduced_size);
284 std::vector<std::pair<NodeIndex, NodeIndex>> closure_arcs;
285 switch (matching_) {
286 case MatchingAlgorithm::MINIMUM_WEIGHT_MATCHING: {
287 auto result = ComputeMinimumWeightMatching(
288 reduced_graph, [this, &reduced_graph,
289 &odd_degree_nodes](CompleteGraph<>::ArcIndex arc) {
290 return costs_(odd_degree_nodes[reduced_graph.Tail(arc)],
291 odd_degree_nodes[reduced_graph.Head(arc)]);
292 });
293 if (!result.ok()) {
294 return false;
295 }
296 result->swap(closure_arcs);
297 break;
298 }
299#if defined(USE_CBC) || defined(USE_SCIP)
300 case MatchingAlgorithm::MINIMUM_WEIGHT_MATCHING_WITH_MIP: {
302 reduced_graph, [this, &reduced_graph,
303 &odd_degree_nodes](CompleteGraph<>::ArcIndex arc) {
304 return costs_(odd_degree_nodes[reduced_graph.Tail(arc)],
305 odd_degree_nodes[reduced_graph.Head(arc)]);
306 });
307 if (!result.ok()) {
308 return false;
309 }
310 result->swap(closure_arcs);
311 break;
312 }
313#endif // defined(USE_CBC) || defined(USE_SCIP)
314 case MatchingAlgorithm::MINIMAL_WEIGHT_MATCHING: {
315 // TODO(user): Cost caching was added and can gain up to 20% but
316 // increases memory usage; see if we can avoid caching.
317 std::vector<ArcIndex> ordered_arcs(reduced_graph.num_arcs());
318 std::vector<CostType> ordered_arc_costs(reduced_graph.num_arcs(), 0);
319 for (const ArcIndex arc : reduced_graph.AllForwardArcs()) {
320 ordered_arcs[arc] = arc;
321 ordered_arc_costs[arc] =
322 costs_(odd_degree_nodes[reduced_graph.Tail(arc)],
323 odd_degree_nodes[reduced_graph.Head(arc)]);
324 }
325 std::sort(ordered_arcs.begin(), ordered_arcs.end(),
326 [&ordered_arc_costs](ArcIndex arc_a, ArcIndex arc_b) {
327 return ordered_arc_costs[arc_a] < ordered_arc_costs[arc_b];
328 });
329 std::vector<bool> touched_nodes(reduced_size, false);
330 for (ArcIndex arc_index = 0; closure_arcs.size() * 2 < reduced_size;
331 ++arc_index) {
332 const ArcIndex arc = ordered_arcs[arc_index];
333 const NodeIndex tail = reduced_graph.Tail(arc);
334 const NodeIndex head = reduced_graph.Head(arc);
335 if (head != tail && !touched_nodes[tail] && !touched_nodes[head]) {
336 touched_nodes[tail] = true;
337 touched_nodes[head] = true;
338 closure_arcs.emplace_back(tail, head);
339 }
340 }
341 break;
342 }
343 }
344 // Build Eulerian path on minimum spanning tree + closing edges from matching
345 // and extract a solution to the Traveling Salesman from the path by skipping
346 // duplicate nodes.
348 num_nodes, closure_arcs.size() + mst.size());
349 for (ArcIndex arc : mst) {
350 egraph.AddArc(graph_.Tail(arc), graph_.Head(arc));
351 }
352 for (const auto arc : closure_arcs) {
353 egraph.AddArc(odd_degree_nodes[arc.first], odd_degree_nodes[arc.second]);
354 }
355 std::vector<bool> touched(num_nodes, false);
356 DCHECK(IsEulerianGraph(egraph));
357 for (const NodeIndex node : BuildEulerianTourFromNode(egraph, 0)) {
358 if (touched[node]) continue;
359 touched[node] = true;
360 tsp_cost_ = SafeAdd(tsp_cost_,
361 tsp_path_.empty() ? 0 : costs_(tsp_path_.back(), node));
362 tsp_path_.push_back(node);
363 }
364 tsp_cost_ =
365 SafeAdd(tsp_cost_, tsp_path_.empty() ? 0 : costs_(tsp_path_.back(), 0));
366 tsp_path_.push_back(0);
367 solved_ = true;
368 return true;
369}
370} // namespace operations_research
371
372#endif // OR_TOOLS_GRAPH_CHRISTOFIDES_H_
std::vector< NodeIndex > TravelingSalesmanPath()
Definition: christofides.h:241
ChristofidesPathSolver(NodeIndex num_nodes, CostFunction costs)
Definition: christofides.h:220
void SetMatchingAlgorithm(MatchingAlgorithm matching)
Definition: christofides.h:61
ArcIndexType AddArc(NodeIndexType tail, NodeIndexType head)
Definition: graph.h:1506
bool IsEulerianGraph(const Graph &graph)
Definition: eulerian_path.h:40
absl::StatusOr< std::vector< std::pair< typename GraphType::NodeIndex, typename GraphType::NodeIndex > > > ComputeMinimumWeightMatchingWithMIP(const GraphType &graph, const WeightFunctionType &weight)
Definition: christofides.h:145
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:109
std::vector< NodeIndex > BuildEulerianTourFromNode(const Graph &graph, NodeIndex root)