OR-Tools  8.2
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"
28#include "ortools/graph/graph.h"
32#include "ortools/linear_solver/linear_solver.pb.h"
34
35namespace operations_research {
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",
193#elif defined(USE_CBC)
194 MPSolver mp_solver("MatchingWithCBC",
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;
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_
#define DCHECK_NE(val1, val2)
Definition: base/logging.h:886
#define DCHECK_GE(val1, val2)
Definition: base/logging.h:889
#define DCHECK(condition)
Definition: base/logging.h:884
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
This mathematical programming (MP) solver class is the main class though which users build and solve ...
void FillSolutionResponseProto(MPSolutionResponse *response) const
Encodes the current solution in a solution response protocol buffer.
ResultStatus
The status of solving the problem.
MPSolverResponseStatus LoadModelFromProto(const MPModelProto &input_model, std::string *error_message)
Loads model from protocol buffer.
ResultStatus Solve()
Solves the problem using the default parameter values.
void AddEdgeWithCost(int tail, int head, int64 cost)
ArcIndexType AddArc(NodeIndexType tail, NodeIndexType head)
Definition: graph.h:1506
SharedResponseManager * response
GRBmodel * model
int int32
int64_t int64
A C++ wrapper that provides a simple and unified interface to several linear programming and mixed in...
CpSolverResponse Solve(const CpModelProto &model_proto)
Solves the given CpModelProto and returns an instance of CpSolverResponse.
The vehicle routing library lets one model and solve generic vehicle routing problems ranging from th...
int64 CapAdd(int64 x, int64 y)
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)
STL namespace.
int64 weight
Definition: pack.cc:509
int64 tail
int64 head