20#ifndef OR_TOOLS_GRAPH_CHRISTOFIDES_H_
21#define OR_TOOLS_GRAPH_CHRISTOFIDES_H_
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"
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"
37using ::util::CompleteGraph;
39template <
typename CostType,
typename ArcIndex = int64,
45 MINIMUM_WEIGHT_MATCHING,
46#if defined(USE_CBC) || defined(USE_SCIP)
47 MINIMUM_WEIGHT_MATCHING_WITH_MIP,
49 MINIMAL_WEIGHT_MATCHING,
82 int64 SafeAdd(int64 a, int64 b) {
90 CompleteGraph<NodeIndex, ArcIndex> graph_;
93 const CostFunction costs_;
99 std::vector<NodeIndex> tsp_path_;
106template <
typename WeightFunctionType,
typename GraphType>
107absl::StatusOr<std::vector<
108 std::pair<typename GraphType::NodeIndex, typename GraphType::NodeIndex>>>
110 const WeightFunctionType& weight) {
113 MinCostPerfectMatching matching(graph.num_nodes());
114 for (
NodeIndex tail : graph.AllNodes()) {
115 for (
const ArcIndex arc : graph.OutgoingArcs(tail)) {
119 matching.AddEdgeWithCost(tail, head, weight(arc));
123 MinCostPerfectMatching::Status status = matching.Solve();
124 if (status != MinCostPerfectMatching::OPTIMAL) {
125 return absl::InvalidArgumentError(
"Perfect matching failed");
127 std::vector<std::pair<NodeIndex, NodeIndex>> match;
128 for (
NodeIndex tail : graph.AllNodes()) {
129 const NodeIndex head = matching.Match(tail);
131 match.emplace_back(tail, head);
137#if defined(USE_CBC) || defined(USE_SCIP)
142template <
typename WeightFunctionType,
typename GraphType>
143absl::StatusOr<std::vector<
144 std::pair<typename GraphType::NodeIndex, typename GraphType::NodeIndex>>>
146 const WeightFunctionType& weight) {
150 model.set_maximize(
false);
155 std::vector<int> variable_indices(graph.num_arcs(), -1);
156 for (
NodeIndex node : graph.AllNodes()) {
158 for (
const ArcIndex arc : graph.OutgoingArcs(node)) {
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));
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);
175 for (
NodeIndex node : graph.AllNodes()) {
176 for (
const ArcIndex arc : graph.OutgoingArcs(node)) {
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);
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);
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");
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));
218 typename CostFunction>
223 costs_(std::move(costs)),
228 typename CostFunction>
230 CostFunction>::TravelingSalesmanCost() {
232 bool const ok = Solve();
239 typename CostFunction>
243 const bool ok = Solve();
250 typename CostFunction>
252 CostFunction>::Solve() {
253 const NodeIndex num_nodes = graph_.num_nodes();
256 if (num_nodes == 1) {
259 if (num_nodes <= 1) {
263 const std::vector<ArcIndex> mst =
265 return costs_(graph_.Tail(arc), graph_.Head(arc));
268 std::vector<NodeIndex> degrees(num_nodes, 0);
270 degrees[graph_.Tail(arc)]++;
271 degrees[graph_.Head(arc)]++;
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);
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;
286 case MatchingAlgorithm::MINIMUM_WEIGHT_MATCHING: {
288 reduced_graph, [
this, &reduced_graph,
290 return costs_(odd_degree_nodes[reduced_graph.Tail(arc)],
291 odd_degree_nodes[reduced_graph.Head(arc)]);
296 result->swap(closure_arcs);
299#if defined(USE_CBC) || defined(USE_SCIP)
300 case MatchingAlgorithm::MINIMUM_WEIGHT_MATCHING_WITH_MIP: {
302 reduced_graph, [
this, &reduced_graph,
304 return costs_(odd_degree_nodes[reduced_graph.Tail(arc)],
305 odd_degree_nodes[reduced_graph.Head(arc)]);
310 result->swap(closure_arcs);
314 case MatchingAlgorithm::MINIMAL_WEIGHT_MATCHING: {
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)]);
325 std::sort(ordered_arcs.begin(), ordered_arcs.end(),
327 return ordered_arc_costs[arc_a] < ordered_arc_costs[arc_b];
329 std::vector<bool> touched_nodes(reduced_size,
false);
330 for (
ArcIndex arc_index = 0; closure_arcs.size() * 2 < reduced_size;
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);
348 num_nodes, closure_arcs.size() + mst.size());
350 egraph.
AddArc(graph_.Tail(arc), graph_.Head(arc));
352 for (
const auto arc : closure_arcs) {
353 egraph.
AddArc(odd_degree_nodes[arc.first], odd_degree_nodes[arc.second]);
355 std::vector<bool> touched(num_nodes,
false);
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);
365 SafeAdd(tsp_cost_, tsp_path_.empty() ? 0 : costs_(tsp_path_.back(), 0));
366 tsp_path_.push_back(0);
std::vector< NodeIndex > TravelingSalesmanPath()
ChristofidesPathSolver(NodeIndex num_nodes, CostFunction costs)
CostType TravelingSalesmanCost()
void SetMatchingAlgorithm(MatchingAlgorithm matching)
ArcIndexType AddArc(NodeIndexType tail, NodeIndexType head)
bool IsEulerianGraph(const Graph &graph)
absl::StatusOr< std::vector< std::pair< typename GraphType::NodeIndex, typename GraphType::NodeIndex > > > ComputeMinimumWeightMatchingWithMIP(const GraphType &graph, const WeightFunctionType &weight)
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)
std::vector< NodeIndex > BuildEulerianTourFromNode(const Graph &graph, NodeIndex root)