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"
26 #include "ortools/base/logging.h"
28 #include "ortools/graph/graph.h"
32 #include "ortools/linear_solver/linear_solver.pb.h"
34 
35 namespace operations_research {
36 
37 using ::util::CompleteGraph;
38 
39 template <typename CostType, typename ArcIndex = int64,
40  typename NodeIndex = int32,
41  typename CostFunction = std::function<CostType(NodeIndex, NodeIndex)>>
43  public:
44  enum class MatchingAlgorithm {
46 #if defined(USE_CBC) || defined(USE_SCIP)
48 #endif // defined(USE_CBC) || defined(USE_SCIP)
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.
106 template <typename WeightFunctionType, typename GraphType>
107 absl::StatusOr<std::vector<
108  std::pair<typename GraphType::NodeIndex, typename GraphType::NodeIndex>>>
109 ComputeMinimumWeightMatching(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.
142 template <typename WeightFunctionType, typename GraphType>
143 absl::StatusOr<std::vector<
144  std::pair<typename GraphType::NodeIndex, typename GraphType::NodeIndex>>>
145 ComputeMinimumWeightMatchingWithMIP(const GraphType& graph,
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 
217 template <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 
227 template <typename CostType, typename ArcIndex, typename NodeIndex,
228  typename CostFunction>
229 CostType 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 
238 template <typename CostType, typename ArcIndex, typename NodeIndex,
239  typename CostFunction>
240 std::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 
249 template <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...
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< typename Graph::ArcIndex > BuildPrimMinimumSpanningTree(const Graph &graph, const ArcValue &arc_value)
absl::StatusOr< std::vector< std::pair< typename GraphType::NodeIndex, typename GraphType::NodeIndex > > > ComputeMinimumWeightMatchingWithMIP(const GraphType &graph, const WeightFunctionType &weight)
Definition: christofides.h:145
int64 CapAdd(int64 x, int64 y)
std::vector< NodeIndex > BuildEulerianTourFromNode(const Graph &graph, NodeIndex root)
bool IsEulerianGraph(const Graph &graph)
Definition: eulerian_path.h:40
int64 weight
Definition: pack.cc:509
int64 tail
int64 head