20 #ifndef OR_TOOLS_GRAPH_CHRISTOFIDES_H_
21 #define OR_TOOLS_GRAPH_CHRISTOFIDES_H_
23 #include "ortools/base/integral_types.h"
24 #include "ortools/base/logging.h"
28 #include "ortools/graph/perfect_matching.h"
29 #include "ortools/linear_solver/linear_solver.h"
30 #include "ortools/linear_solver/linear_solver.pb.h"
31 #include "ortools/util/saturated_arithmetic.h"
35 using ::util::CompleteGraph;
37 template <
typename CostType,
typename ArcIndex = int64,
44 #if defined(USE_CBC) || defined(USE_SCIP)
46 #endif // defined(USE_CBC) || defined(USE_SCIP)
79 int64 SafeAdd(int64 a, int64 b) {
87 CompleteGraph<NodeIndex, ArcIndex> graph_;
90 const CostFunction costs_;
96 std::vector<NodeIndex> tsp_path_;
103 template <
typename WeightFunctionType,
typename GraphType>
105 std::pair<typename GraphType::NodeIndex, typename GraphType::NodeIndex>>
107 const WeightFunctionType& weight) {
110 MinCostPerfectMatching matching(graph.num_nodes());
111 for (
NodeIndex tail : graph.AllNodes()) {
112 for (
const ArcIndex arc : graph.OutgoingArcs(tail)) {
116 matching.AddEdgeWithCost(tail, head, weight(arc));
120 MinCostPerfectMatching::Status status = matching.Solve();
121 DCHECK_EQ(status, MinCostPerfectMatching::OPTIMAL);
122 std::vector<std::pair<NodeIndex, NodeIndex>> match;
123 for (
NodeIndex tail : graph.AllNodes()) {
124 const NodeIndex head = matching.Match(tail);
126 match.emplace_back(tail, head);
132 #if defined(USE_CBC) || defined(USE_SCIP)
137 template <
typename WeightFunctionType,
typename GraphType>
139 std::pair<typename GraphType::NodeIndex, typename GraphType::NodeIndex>>
141 const WeightFunctionType& weight) {
145 model.set_maximize(
false);
150 std::vector<int> variable_indices(graph.num_arcs(), -1);
151 for (
NodeIndex node : graph.AllNodes()) {
153 for (
const ArcIndex arc : graph.OutgoingArcs(node)) {
156 variable_indices[arc] = model.variable_size();
157 MPVariableProto*
const arc_var = model.add_variable();
158 arc_var->set_lower_bound(0);
159 arc_var->set_upper_bound(1);
160 arc_var->set_is_integer(
true);
161 arc_var->set_objective_coefficient(weight(arc));
166 MPConstraintProto*
const one_of_ct = model.add_constraint();
167 one_of_ct->set_lower_bound(1);
168 one_of_ct->set_upper_bound(1);
170 for (
NodeIndex node : graph.AllNodes()) {
171 for (
const ArcIndex arc : graph.OutgoingArcs(node)) {
174 const int arc_var = variable_indices[arc];
175 DCHECK_GE(arc_var, 0);
176 MPConstraintProto* one_of_ct = model.mutable_constraint(node);
177 one_of_ct->add_var_index(arc_var);
178 one_of_ct->add_coefficient(1);
179 one_of_ct = model.mutable_constraint(head);
180 one_of_ct->add_var_index(arc_var);
181 one_of_ct->add_coefficient(1);
185 #if defined(USE_SCIP)
186 MPSolver mp_solver(
"MatchingWithSCIP",
187 MPSolver::SCIP_MIXED_INTEGER_PROGRAMMING);
188 #elif defined(USE_CBC)
189 MPSolver mp_solver(
"MatchingWithCBC",
190 MPSolver::CBC_MIXED_INTEGER_PROGRAMMING);
193 mp_solver.LoadModelFromProto(model, &error);
194 MPSolver::ResultStatus status = mp_solver.Solve();
195 CHECK_EQ(status, MPSolver::OPTIMAL);
196 MPSolutionResponse response;
197 mp_solver.FillSolutionResponseProto(&response);
198 std::vector<std::pair<NodeIndex, NodeIndex>> matching;
199 for (
ArcIndex arc = 0; arc < variable_indices.size(); ++arc) {
200 const int arc_var = variable_indices[arc];
201 if (arc_var >= 0 && response.variable_value(arc_var) > .9) {
202 DCHECK_GE(response.variable_value(arc_var), 1.0 - 1e-4);
203 matching.emplace_back(graph.Tail(arc), graph.Head(arc));
208 #endif // defined(USE_CBC) || defined(USE_SCIP)
211 typename CostFunction>
216 costs_(std::move(costs)),
221 typename CostFunction>
223 CostFunction>::TravelingSalesmanCost() {
231 typename CostFunction>
241 typename CostFunction>
243 CostFunction>::Solve() {
244 const NodeIndex num_nodes = graph_.num_nodes();
247 if (num_nodes == 1) {
250 if (num_nodes <= 1) {
254 const std::vector<ArcIndex> mst =
256 return costs_(graph_.Tail(arc), graph_.Head(arc));
259 std::vector<NodeIndex> degrees(num_nodes, 0);
261 degrees[graph_.Tail(arc)]++;
262 degrees[graph_.Head(arc)]++;
264 std::vector<NodeIndex> odd_degree_nodes;
265 for (
int i = 0; i < degrees.size(); ++i) {
266 if (degrees[i] % 2 != 0) {
267 odd_degree_nodes.push_back(i);
272 const NodeIndex reduced_size = odd_degree_nodes.size();
273 DCHECK_NE(0, reduced_size);
274 CompleteGraph<NodeIndex, ArcIndex> reduced_graph(reduced_size);
275 std::vector<std::pair<NodeIndex, NodeIndex>> closure_arcs;
277 case MatchingAlgorithm::MINIMUM_WEIGHT_MATCHING: {
279 reduced_graph, [
this, &reduced_graph,
281 return costs_(odd_degree_nodes[reduced_graph.Tail(arc)],
282 odd_degree_nodes[reduced_graph.Head(arc)]);
286 #if defined(USE_CBC) || defined(USE_SCIP)
287 case MatchingAlgorithm::MINIMUM_WEIGHT_MATCHING_WITH_MIP: {
289 reduced_graph, [
this, &reduced_graph,
291 return costs_(odd_degree_nodes[reduced_graph.Tail(arc)],
292 odd_degree_nodes[reduced_graph.Head(arc)]);
296 #endif // defined(USE_CBC) || defined(USE_SCIP)
297 case MatchingAlgorithm::MINIMAL_WEIGHT_MATCHING: {
300 std::vector<ArcIndex> ordered_arcs(reduced_graph.num_arcs());
301 std::vector<CostType> ordered_arc_costs(reduced_graph.num_arcs(), 0);
302 for (
const ArcIndex arc : reduced_graph.AllForwardArcs()) {
303 ordered_arcs[arc] = arc;
304 ordered_arc_costs[arc] =
305 costs_(odd_degree_nodes[reduced_graph.Tail(arc)],
306 odd_degree_nodes[reduced_graph.Head(arc)]);
308 std::sort(ordered_arcs.begin(), ordered_arcs.end(),
310 return ordered_arc_costs[arc_a] < ordered_arc_costs[arc_b];
312 std::vector<bool> touched_nodes(reduced_size,
false);
313 for (
ArcIndex arc_index = 0; closure_arcs.size() * 2 < reduced_size;
315 const ArcIndex arc = ordered_arcs[arc_index];
316 const NodeIndex tail = reduced_graph.Tail(arc);
317 const NodeIndex head = reduced_graph.Head(arc);
318 if (head != tail && !touched_nodes[tail] && !touched_nodes[head]) {
319 touched_nodes[tail] =
true;
320 touched_nodes[head] =
true;
321 closure_arcs.emplace_back(tail, head);
331 num_nodes, closure_arcs.size() + mst.size());
333 egraph.AddArc(graph_.Tail(arc), graph_.Head(arc));
335 for (
const auto arc : closure_arcs) {
336 egraph.AddArc(odd_degree_nodes[arc.first], odd_degree_nodes[arc.second]);
338 std::vector<bool> touched(num_nodes,
false);
341 if (touched[node])
continue;
342 touched[node] =
true;
343 tsp_cost_ = SafeAdd(tsp_cost_,
344 tsp_path_.empty() ? 0 : costs_(tsp_path_.back(), node));
345 tsp_path_.push_back(node);
348 SafeAdd(tsp_cost_, tsp_path_.empty() ? 0 : costs_(tsp_path_.back(), 0));
349 tsp_path_.push_back(0);
354 #endif // OR_TOOLS_GRAPH_CHRISTOFIDES_H_