Shortcuts on this page
r m x toggle line displays
j k next/prev highlighted chunk
0 (zero) top of page
1 (one) first highlighted chunk
1"""
2=================================
3Travelling Salesman Problem (TSP)
4=================================
6Implementation of approximate algorithms
7for solving and approximating the TSP problem.
9Categories of algorithms which are implemented:
11- Christofides (provides a 3/2-approximation of TSP)
12- Greedy
13- Simulated Annealing (SA)
14- Threshold Accepting (TA)
15- Asadpour Asymmetric Traveling Salesman Algorithm
17The Travelling Salesman Problem tries to find, given the weight
18(distance) between all points where a salesman has to visit, the
19route so that:
21- The total distance (cost) which the salesman travels is minimized.
22- The salesman returns to the starting point.
23- Note that for a complete graph, the salesman visits each point once.
25The function `travelling_salesman_problem` allows for incomplete
26graphs by finding all-pairs shortest paths, effectively converting
27the problem to a complete graph problem. It calls one of the
28approximate methods on that problem and then converts the result
29back to the original graph using the previously found shortest paths.
31TSP is an NP-hard problem in combinatorial optimization,
32important in operations research and theoretical computer science.
34http://en.wikipedia.org/wiki/Travelling_salesman_problem
35"""
37import math
39import networkx as nx
40from networkx.algorithms.tree.mst import random_spanning_tree
41from networkx.utils import not_implemented_for, pairwise, py_random_state
43__all__ = [
44 "traveling_salesman_problem",
45 "christofides",
46 "asadpour_atsp",
47 "greedy_tsp",
48 "simulated_annealing_tsp",
49 "threshold_accepting_tsp",
50]
53def swap_two_nodes(soln, seed):
54 """Swap two nodes in `soln` to give a neighbor solution.
56 Parameters
57 ----------
58 soln : list of nodes
59 Current cycle of nodes
61 seed : integer, random_state, or None (default)
62 Indicator of random number generation state.
63 See :ref:`Randomness<randomness>`.
65 Returns
66 -------
67 list
68 The solution after move is applied. (A neighbor solution.)
70 Notes
71 -----
72 This function assumes that the incoming list `soln` is a cycle
73 (that the first and last element are the same) and also that
74 we don't want any move to change the first node in the list
75 (and thus not the last node either).
77 The input list is changed as well as returned. Make a copy if needed.
79 See Also
80 --------
81 move_one_node
82 """
83 a, b = seed.sample(range(1, len(soln) - 1), k=2)
84 soln[a], soln[b] = soln[b], soln[a]
85 return soln
88def move_one_node(soln, seed):
89 """Move one node to another position to give a neighbor solution.
91 The node to move and the position to move to are chosen randomly.
92 The first and last nodes are left untouched as soln must be a cycle
93 starting at that node.
95 Parameters
96 ----------
97 soln : list of nodes
98 Current cycle of nodes
100 seed : integer, random_state, or None (default)
101 Indicator of random number generation state.
102 See :ref:`Randomness<randomness>`.
104 Returns
105 -------
106 list
107 The solution after move is applied. (A neighbor solution.)
109 Notes
110 -----
111 This function assumes that the incoming list `soln` is a cycle
112 (that the first and last element are the same) and also that
113 we don't want any move to change the first node in the list
114 (and thus not the last node either).
116 The input list is changed as well as returned. Make a copy if needed.
118 See Also
119 --------
120 swap_two_nodes
121 """
122 a, b = seed.sample(range(1, len(soln) - 1), k=2)
123 soln.insert(b, soln.pop(a))
124 return soln
127@not_implemented_for("directed")
128@nx._dispatchable(edge_attrs="weight")
129def christofides(G, weight="weight", tree=None):
130 """Approximate a solution of the traveling salesman problem
132 Compute a 3/2-approximation of the traveling salesman problem
133 in a complete undirected graph using Christofides [1]_ algorithm.
135 Parameters
136 ----------
137 G : Graph
138 `G` should be a complete weighted undirected graph.
139 The distance between all pairs of nodes should be included.
141 weight : string, optional (default="weight")
142 Edge data key corresponding to the edge weight.
143 If any edge does not have this attribute the weight is set to 1.
145 tree : NetworkX graph or None (default: None)
146 A minimum spanning tree of G. Or, if None, the minimum spanning
147 tree is computed using :func:`networkx.minimum_spanning_tree`
149 Returns
150 -------
151 list
152 List of nodes in `G` along a cycle with a 3/2-approximation of
153 the minimal Hamiltonian cycle.
155 References
156 ----------
157 .. [1] Christofides, Nicos. "Worst-case analysis of a new heuristic for
158 the travelling salesman problem." No. RR-388. Carnegie-Mellon Univ
159 Pittsburgh Pa Management Sciences Research Group, 1976.
160 """
161 # Remove selfloops if necessary
162 loop_nodes = nx.nodes_with_selfloops(G)
163 try:
164 node = next(loop_nodes)
165 except StopIteration:
166 pass
167 else:
168 G = G.copy()
169 G.remove_edge(node, node)
170 G.remove_edges_from((n, n) for n in loop_nodes)
171 # Check that G is a complete graph
172 N = len(G) - 1
173 # This check ignores selfloops which is what we want here.
174 if any(len(nbrdict) != N for n, nbrdict in G.adj.items()):
175 raise nx.NetworkXError("G must be a complete graph.")
177 if tree is None:
178 tree = nx.minimum_spanning_tree(G, weight=weight)
179 L = G.copy()
180 L.remove_nodes_from([v for v, degree in tree.degree if not (degree % 2)])
181 MG = nx.MultiGraph()
182 MG.add_edges_from(tree.edges)
183 edges = nx.min_weight_matching(L, weight=weight)
184 MG.add_edges_from(edges)
185 return _shortcutting(nx.eulerian_circuit(MG))
188def _shortcutting(circuit):
189 """Remove duplicate nodes in the path"""
190 nodes = []
191 for u, v in circuit:
192 if v in nodes:
193 continue
194 if not nodes:
195 nodes.append(u)
196 nodes.append(v)
197 nodes.append(nodes[0])
198 return nodes
201@nx._dispatchable(edge_attrs="weight")
202def traveling_salesman_problem(
203 G, weight="weight", nodes=None, cycle=True, method=None, **kwargs
204):
205 """Find the shortest path in `G` connecting specified nodes
207 This function allows approximate solution to the traveling salesman
208 problem on networks that are not complete graphs and/or where the
209 salesman does not need to visit all nodes.
211 This function proceeds in two steps. First, it creates a complete
212 graph using the all-pairs shortest_paths between nodes in `nodes`.
213 Edge weights in the new graph are the lengths of the paths
214 between each pair of nodes in the original graph.
215 Second, an algorithm (default: `christofides` for undirected and
216 `asadpour_atsp` for directed) is used to approximate the minimal Hamiltonian
217 cycle on this new graph. The available algorithms are:
219 - christofides
220 - greedy_tsp
221 - simulated_annealing_tsp
222 - threshold_accepting_tsp
223 - asadpour_atsp
225 Once the Hamiltonian Cycle is found, this function post-processes to
226 accommodate the structure of the original graph. If `cycle` is ``False``,
227 the biggest weight edge is removed to make a Hamiltonian path.
228 Then each edge on the new complete graph used for that analysis is
229 replaced by the shortest_path between those nodes on the original graph.
230 If the input graph `G` includes edges with weights that do not adhere to
231 the triangle inequality, such as when `G` is not a complete graph (i.e
232 length of non-existent edges is infinity), then the returned path may
233 contain some repeating nodes (other than the starting node).
235 Parameters
236 ----------
237 G : NetworkX graph
238 A possibly weighted graph
240 nodes : collection of nodes (default=G.nodes)
241 collection (list, set, etc.) of nodes to visit
243 weight : string, optional (default="weight")
244 Edge data key corresponding to the edge weight.
245 If any edge does not have this attribute the weight is set to 1.
247 cycle : bool (default: True)
248 Indicates whether a cycle should be returned, or a path.
249 Note: the cycle is the approximate minimal cycle.
250 The path simply removes the biggest edge in that cycle.
252 method : function (default: None)
253 A function that returns a cycle on all nodes and approximates
254 the solution to the traveling salesman problem on a complete
255 graph. The returned cycle is then used to find a corresponding
256 solution on `G`. `method` should be callable; take inputs
257 `G`, and `weight`; and return a list of nodes along the cycle.
259 Provided options include :func:`christofides`, :func:`greedy_tsp`,
260 :func:`simulated_annealing_tsp` and :func:`threshold_accepting_tsp`.
262 If `method is None`: use :func:`christofides` for undirected `G` and
263 :func:`asadpour_atsp` for directed `G`.
265 **kwargs : dict
266 Other keyword arguments to be passed to the `method` function passed in.
268 Returns
269 -------
270 list
271 List of nodes in `G` along a path with an approximation of the minimal
272 path through `nodes`.
274 Raises
275 ------
276 NetworkXError
277 If `G` is a directed graph it has to be strongly connected or the
278 complete version cannot be generated.
280 Examples
281 --------
282 >>> tsp = nx.approximation.traveling_salesman_problem
283 >>> G = nx.cycle_graph(9)
284 >>> G[4][5]["weight"] = 5 # all other weights are 1
285 >>> tsp(G, nodes=[3, 6])
286 [3, 2, 1, 0, 8, 7, 6, 7, 8, 0, 1, 2, 3]
287 >>> path = tsp(G, cycle=False)
288 >>> path in ([4, 3, 2, 1, 0, 8, 7, 6, 5], [5, 6, 7, 8, 0, 1, 2, 3, 4])
289 True
291 While no longer required, you can still build (curry) your own function
292 to provide parameter values to the methods.
294 >>> SA_tsp = nx.approximation.simulated_annealing_tsp
295 >>> method = lambda G, weight: SA_tsp(G, "greedy", weight=weight, temp=500)
296 >>> path = tsp(G, cycle=False, method=method)
297 >>> path in ([4, 3, 2, 1, 0, 8, 7, 6, 5], [5, 6, 7, 8, 0, 1, 2, 3, 4])
298 True
300 Otherwise, pass other keyword arguments directly into the tsp function.
302 >>> path = tsp(
303 ... G,
304 ... cycle=False,
305 ... method=nx.approximation.simulated_annealing_tsp,
306 ... init_cycle="greedy",
307 ... temp=500,
308 ... )
309 >>> path in ([4, 3, 2, 1, 0, 8, 7, 6, 5], [5, 6, 7, 8, 0, 1, 2, 3, 4])
310 True
311 """
312 if method is None:
313 if G.is_directed():
314 method = asadpour_atsp
315 else:
316 method = christofides
317 if nodes is None:
318 nodes = list(G.nodes)
320 dist = {}
321 path = {}
322 for n, (d, p) in nx.all_pairs_dijkstra(G, weight=weight):
323 dist[n] = d
324 path[n] = p
326 if G.is_directed():
327 # If the graph is not strongly connected, raise an exception
328 if not nx.is_strongly_connected(G):
329 raise nx.NetworkXError("G is not strongly connected")
330 GG = nx.DiGraph()
331 else:
332 GG = nx.Graph()
333 for u in nodes:
334 for v in nodes:
335 if u == v:
336 continue
337 # Ensure that the weight attribute on GG has the
338 # same name as the input graph
339 GG.add_edge(u, v, **{weight: dist[u][v]})
341 best_GG = method(GG, weight=weight, **kwargs)
343 if not cycle:
344 # find and remove the biggest edge
345 (u, v) = max(pairwise(best_GG), key=lambda x: dist[x[0]][x[1]])
346 pos = best_GG.index(u) + 1
347 while best_GG[pos] != v:
348 pos = best_GG[pos:].index(u) + 1
349 best_GG = best_GG[pos:-1] + best_GG[:pos]
351 best_path = []
352 for u, v in pairwise(best_GG):
353 best_path.extend(path[u][v][:-1])
354 best_path.append(v)
355 return best_path
358@not_implemented_for("undirected")
359@py_random_state(2)
360@nx._dispatchable(edge_attrs="weight", mutates_input=True)
361def asadpour_atsp(G, weight="weight", seed=None, source=None):
362 """
363 Returns an approximate solution to the traveling salesman problem.
365 This approximate solution is one of the best known approximations for the
366 asymmetric traveling salesman problem developed by Asadpour et al,
367 [1]_. The algorithm first solves the Held-Karp relaxation to find a lower
368 bound for the weight of the cycle. Next, it constructs an exponential
369 distribution of undirected spanning trees where the probability of an
370 edge being in the tree corresponds to the weight of that edge using a
371 maximum entropy rounding scheme. Next we sample that distribution
372 $2 \\lceil \\ln n \\rceil$ times and save the minimum sampled tree once the
373 direction of the arcs is added back to the edges. Finally, we augment
374 then short circuit that graph to find the approximate tour for the
375 salesman.
377 Parameters
378 ----------
379 G : nx.DiGraph
380 The graph should be a complete weighted directed graph. The
381 distance between all paris of nodes should be included and the triangle
382 inequality should hold. That is, the direct edge between any two nodes
383 should be the path of least cost.
385 weight : string, optional (default="weight")
386 Edge data key corresponding to the edge weight.
387 If any edge does not have this attribute the weight is set to 1.
389 seed : integer, random_state, or None (default)
390 Indicator of random number generation state.
391 See :ref:`Randomness<randomness>`.
393 source : node label (default=`None`)
394 If given, return the cycle starting and ending at the given node.
396 Returns
397 -------
398 cycle : list of nodes
399 Returns the cycle (list of nodes) that a salesman can follow to minimize
400 the total weight of the trip.
402 Raises
403 ------
404 NetworkXError
405 If `G` is not complete or has less than two nodes, the algorithm raises
406 an exception.
408 NetworkXError
409 If `source` is not `None` and is not a node in `G`, the algorithm raises
410 an exception.
412 NetworkXNotImplemented
413 If `G` is an undirected graph.
415 References
416 ----------
417 .. [1] A. Asadpour, M. X. Goemans, A. Madry, S. O. Gharan, and A. Saberi,
418 An o(log n/log log n)-approximation algorithm for the asymmetric
419 traveling salesman problem, Operations research, 65 (2017),
420 pp. 1043–1061
422 Examples
423 --------
424 >>> import networkx as nx
425 >>> import networkx.algorithms.approximation as approx
426 >>> G = nx.complete_graph(3, create_using=nx.DiGraph)
427 >>> nx.set_edge_attributes(
428 ... G,
429 ... {(0, 1): 2, (1, 2): 2, (2, 0): 2, (0, 2): 1, (2, 1): 1, (1, 0): 1},
430 ... "weight",
431 ... )
432 >>> tour = approx.asadpour_atsp(G, source=0)
433 >>> tour
434 [0, 2, 1, 0]
435 """
436 from math import ceil, exp
437 from math import log as ln
439 # Check that G is a complete graph
440 N = len(G) - 1
441 if N < 1:
442 raise nx.NetworkXError("G must have at least two nodes")
443 # This check ignores selfloops which is what we want here.
444 if any(len(nbrdict) - (n in nbrdict) != N for n, nbrdict in G.adj.items()):
445 raise nx.NetworkXError("G is not a complete DiGraph")
446 # Check that the source vertex, if given, is in the graph
447 if source is not None and source not in G.nodes:
448 raise nx.NetworkXError("Given source node not in G.")
449 # handle 2 node case
450 if N == 1:
451 if source is None:
452 return list(G)
453 return [source, next(n for n in G if n != source)]
455 opt_hk, z_star = held_karp_ascent(G, weight)
457 # Test to see if the ascent method found an integer solution or a fractional
458 # solution. If it is integral then z_star is a nx.Graph, otherwise it is
459 # a dict
460 if not isinstance(z_star, dict):
461 # Here we are using the shortcutting method to go from the list of edges
462 # returned from eulerian_circuit to a list of nodes
463 return _shortcutting(nx.eulerian_circuit(z_star, source=source))
465 # Create the undirected support of z_star
466 z_support = nx.MultiGraph()
467 for u, v in z_star:
468 if (u, v) not in z_support.edges:
469 edge_weight = min(G[u][v][weight], G[v][u][weight])
470 z_support.add_edge(u, v, **{weight: edge_weight})
472 # Create the exponential distribution of spanning trees
473 gamma = spanning_tree_distribution(z_support, z_star)
475 # Write the lambda values to the edges of z_support
476 z_support = nx.Graph(z_support)
477 lambda_dict = {(u, v): exp(gamma[(u, v)]) for u, v in z_support.edges()}
478 nx.set_edge_attributes(z_support, lambda_dict, "weight")
479 del gamma, lambda_dict
481 # Sample 2 * ceil( ln(n) ) spanning trees and record the minimum one
482 minimum_sampled_tree = None
483 minimum_sampled_tree_weight = math.inf
484 for _ in range(2 * ceil(ln(G.number_of_nodes()))):
485 sampled_tree = random_spanning_tree(z_support, "weight", seed=seed)
486 sampled_tree_weight = sampled_tree.size(weight)
487 if sampled_tree_weight < minimum_sampled_tree_weight:
488 minimum_sampled_tree = sampled_tree.copy()
489 minimum_sampled_tree_weight = sampled_tree_weight
491 # Orient the edges in that tree to keep the cost of the tree the same.
492 t_star = nx.MultiDiGraph()
493 for u, v, d in minimum_sampled_tree.edges(data=weight):
494 if d == G[u][v][weight]:
495 t_star.add_edge(u, v, **{weight: d})
496 else:
497 t_star.add_edge(v, u, **{weight: d})
499 # Find the node demands needed to neutralize the flow of t_star in G
500 node_demands = {n: t_star.out_degree(n) - t_star.in_degree(n) for n in t_star}
501 nx.set_node_attributes(G, node_demands, "demand")
503 # Find the min_cost_flow
504 flow_dict = nx.min_cost_flow(G, "demand")
506 # Build the flow into t_star
507 for source, values in flow_dict.items():
508 for target in values:
509 if (source, target) not in t_star.edges and values[target] > 0:
510 # IF values[target] > 0 we have to add that many edges
511 for _ in range(values[target]):
512 t_star.add_edge(source, target)
514 # Return the shortcut eulerian circuit
515 circuit = nx.eulerian_circuit(t_star, source=source)
516 return _shortcutting(circuit)
519@nx._dispatchable(edge_attrs="weight", mutates_input=True, returns_graph=True)
520def held_karp_ascent(G, weight="weight"):
521 """
522 Minimizes the Held-Karp relaxation of the TSP for `G`
524 Solves the Held-Karp relaxation of the input complete digraph and scales
525 the output solution for use in the Asadpour [1]_ ASTP algorithm.
527 The Held-Karp relaxation defines the lower bound for solutions to the
528 ATSP, although it does return a fractional solution. This is used in the
529 Asadpour algorithm as an initial solution which is later rounded to a
530 integral tree within the spanning tree polytopes. This function solves
531 the relaxation with the branch and bound method in [2]_.
533 Parameters
534 ----------
535 G : nx.DiGraph
536 The graph should be a complete weighted directed graph.
537 The distance between all paris of nodes should be included.
539 weight : string, optional (default="weight")
540 Edge data key corresponding to the edge weight.
541 If any edge does not have this attribute the weight is set to 1.
543 Returns
544 -------
545 OPT : float
546 The cost for the optimal solution to the Held-Karp relaxation
547 z : dict or nx.Graph
548 A symmetrized and scaled version of the optimal solution to the
549 Held-Karp relaxation for use in the Asadpour algorithm.
551 If an integral solution is found, then that is an optimal solution for
552 the ATSP problem and that is returned instead.
554 References
555 ----------
556 .. [1] A. Asadpour, M. X. Goemans, A. Madry, S. O. Gharan, and A. Saberi,
557 An o(log n/log log n)-approximation algorithm for the asymmetric
558 traveling salesman problem, Operations research, 65 (2017),
559 pp. 1043–1061
561 .. [2] M. Held, R. M. Karp, The traveling-salesman problem and minimum
562 spanning trees, Operations Research, 1970-11-01, Vol. 18 (6),
563 pp.1138-1162
564 """
565 import numpy as np
566 import scipy as sp
568 def k_pi():
569 """
570 Find the set of minimum 1-Arborescences for G at point pi.
572 Returns
573 -------
574 Set
575 The set of minimum 1-Arborescences
576 """
577 # Create a copy of G without vertex 1.
578 G_1 = G.copy()
579 minimum_1_arborescences = set()
580 minimum_1_arborescence_weight = math.inf
582 # node is node '1' in the Held and Karp paper
583 n = next(G.__iter__())
584 G_1.remove_node(n)
586 # Iterate over the spanning arborescences of the graph until we know
587 # that we have found the minimum 1-arborescences. My proposed strategy
588 # is to find the most extensive root to connect to from 'node 1' and
589 # the least expensive one. We then iterate over arborescences until
590 # the cost of the basic arborescence is the cost of the minimum one
591 # plus the difference between the most and least expensive roots,
592 # that way the cost of connecting 'node 1' will by definition not by
593 # minimum
594 min_root = {"node": None, weight: math.inf}
595 max_root = {"node": None, weight: -math.inf}
596 for u, v, d in G.edges(n, data=True):
597 if d[weight] < min_root[weight]:
598 min_root = {"node": v, weight: d[weight]}
599 if d[weight] > max_root[weight]:
600 max_root = {"node": v, weight: d[weight]}
602 min_in_edge = min(G.in_edges(n, data=True), key=lambda x: x[2][weight])
603 min_root[weight] = min_root[weight] + min_in_edge[2][weight]
604 max_root[weight] = max_root[weight] + min_in_edge[2][weight]
606 min_arb_weight = math.inf
607 for arb in nx.ArborescenceIterator(G_1):
608 arb_weight = arb.size(weight)
609 if min_arb_weight == math.inf:
610 min_arb_weight = arb_weight
611 elif arb_weight > min_arb_weight + max_root[weight] - min_root[weight]:
612 break
613 # We have to pick the root node of the arborescence for the out
614 # edge of the first vertex as that is the only node without an
615 # edge directed into it.
616 for N, deg in arb.in_degree:
617 if deg == 0:
618 # root found
619 arb.add_edge(n, N, **{weight: G[n][N][weight]})
620 arb_weight += G[n][N][weight]
621 break
623 # We can pick the minimum weight in-edge for the vertex with
624 # a cycle. If there are multiple edges with the same, minimum
625 # weight, We need to add all of them.
626 #
627 # Delete the edge (N, v) so that we cannot pick it.
628 edge_data = G[N][n]
629 G.remove_edge(N, n)
630 min_weight = min(G.in_edges(n, data=weight), key=lambda x: x[2])[2]
631 min_edges = [
632 (u, v, d) for u, v, d in G.in_edges(n, data=weight) if d == min_weight
633 ]
634 for u, v, d in min_edges:
635 new_arb = arb.copy()
636 new_arb.add_edge(u, v, **{weight: d})
637 new_arb_weight = arb_weight + d
638 # Check to see the weight of the arborescence, if it is a
639 # new minimum, clear all of the old potential minimum
640 # 1-arborescences and add this is the only one. If its
641 # weight is above the known minimum, do not add it.
642 if new_arb_weight < minimum_1_arborescence_weight:
643 minimum_1_arborescences.clear()
644 minimum_1_arborescence_weight = new_arb_weight
645 # We have a 1-arborescence, add it to the set
646 if new_arb_weight == minimum_1_arborescence_weight:
647 minimum_1_arborescences.add(new_arb)
648 G.add_edge(N, n, **edge_data)
650 return minimum_1_arborescences
652 def direction_of_ascent():
653 """
654 Find the direction of ascent at point pi.
656 See [1]_ for more information.
658 Returns
659 -------
660 dict
661 A mapping from the nodes of the graph which represents the direction
662 of ascent.
664 References
665 ----------
666 .. [1] M. Held, R. M. Karp, The traveling-salesman problem and minimum
667 spanning trees, Operations Research, 1970-11-01, Vol. 18 (6),
668 pp.1138-1162
669 """
670 # 1. Set d equal to the zero n-vector.
671 d = {}
672 for n in G:
673 d[n] = 0
674 del n
675 # 2. Find a 1-Arborescence T^k such that k is in K(pi, d).
676 minimum_1_arborescences = k_pi()
677 while True:
678 # Reduce K(pi) to K(pi, d)
679 # Find the arborescence in K(pi) which increases the lest in
680 # direction d
681 min_k_d_weight = math.inf
682 min_k_d = None
683 for arborescence in minimum_1_arborescences:
684 weighted_cost = 0
685 for n, deg in arborescence.degree:
686 weighted_cost += d[n] * (deg - 2)
687 if weighted_cost < min_k_d_weight:
688 min_k_d_weight = weighted_cost
689 min_k_d = arborescence
691 # 3. If sum of d_i * v_{i, k} is greater than zero, terminate
692 if min_k_d_weight > 0:
693 return d, min_k_d
694 # 4. d_i = d_i + v_{i, k}
695 for n, deg in min_k_d.degree:
696 d[n] += deg - 2
697 # Check that we do not need to terminate because the direction
698 # of ascent does not exist. This is done with linear
699 # programming.
700 c = np.full(len(minimum_1_arborescences), -1, dtype=int)
701 a_eq = np.empty((len(G) + 1, len(minimum_1_arborescences)), dtype=int)
702 b_eq = np.zeros(len(G) + 1, dtype=int)
703 b_eq[len(G)] = 1
704 for arb_count, arborescence in enumerate(minimum_1_arborescences):
705 n_count = len(G) - 1
706 for n, deg in arborescence.degree:
707 a_eq[n_count][arb_count] = deg - 2
708 n_count -= 1
709 a_eq[len(G)][arb_count] = 1
710 program_result = sp.optimize.linprog(
711 c, A_eq=a_eq, b_eq=b_eq, method="highs-ipm"
712 )
713 # If the constants exist, then the direction of ascent doesn't
714 if program_result.success:
715 # There is no direction of ascent
716 return None, minimum_1_arborescences
718 # 5. GO TO 2
720 def find_epsilon(k, d):
721 """
722 Given the direction of ascent at pi, find the maximum distance we can go
723 in that direction.
725 Parameters
726 ----------
727 k_xy : set
728 The set of 1-arborescences which have the minimum rate of increase
729 in the direction of ascent
731 d : dict
732 The direction of ascent
734 Returns
735 -------
736 float
737 The distance we can travel in direction `d`
738 """
739 min_epsilon = math.inf
740 for e_u, e_v, e_w in G.edges(data=weight):
741 if (e_u, e_v) in k.edges:
742 continue
743 # Now, I have found a condition which MUST be true for the edges to
744 # be a valid substitute. The edge in the graph which is the
745 # substitute is the one with the same terminated end. This can be
746 # checked rather simply.
747 #
748 # Find the edge within k which is the substitute. Because k is a
749 # 1-arborescence, we know that they is only one such edges
750 # leading into every vertex.
751 if len(k.in_edges(e_v, data=weight)) > 1:
752 raise Exception
753 sub_u, sub_v, sub_w = next(k.in_edges(e_v, data=weight).__iter__())
754 k.add_edge(e_u, e_v, **{weight: e_w})
755 k.remove_edge(sub_u, sub_v)
756 if (
757 max(d for n, d in k.in_degree()) <= 1
758 and len(G) == k.number_of_edges()
759 and nx.is_weakly_connected(k)
760 ):
761 # Ascent method calculation
762 if d[sub_u] == d[e_u] or sub_w == e_w:
763 # Revert to the original graph
764 k.remove_edge(e_u, e_v)
765 k.add_edge(sub_u, sub_v, **{weight: sub_w})
766 continue
767 epsilon = (sub_w - e_w) / (d[e_u] - d[sub_u])
768 if 0 < epsilon < min_epsilon:
769 min_epsilon = epsilon
770 # Revert to the original graph
771 k.remove_edge(e_u, e_v)
772 k.add_edge(sub_u, sub_v, **{weight: sub_w})
774 return min_epsilon
776 # I have to know that the elements in pi correspond to the correct elements
777 # in the direction of ascent, even if the node labels are not integers.
778 # Thus, I will use dictionaries to made that mapping.
779 pi_dict = {}
780 for n in G:
781 pi_dict[n] = 0
782 del n
783 original_edge_weights = {}
784 for u, v, d in G.edges(data=True):
785 original_edge_weights[(u, v)] = d[weight]
786 dir_ascent, k_d = direction_of_ascent()
787 while dir_ascent is not None:
788 max_distance = find_epsilon(k_d, dir_ascent)
789 for n, v in dir_ascent.items():
790 pi_dict[n] += max_distance * v
791 for u, v, d in G.edges(data=True):
792 d[weight] = original_edge_weights[(u, v)] + pi_dict[u]
793 dir_ascent, k_d = direction_of_ascent()
794 nx._clear_cache(G)
795 # k_d is no longer an individual 1-arborescence but rather a set of
796 # minimal 1-arborescences at the maximum point of the polytope and should
797 # be reflected as such
798 k_max = k_d
800 # Search for a cycle within k_max. If a cycle exists, return it as the
801 # solution
802 for k in k_max:
803 if len([n for n in k if k.degree(n) == 2]) == G.order():
804 # Tour found
805 # TODO: this branch does not restore original_edge_weights of G!
806 return k.size(weight), k
808 # Write the original edge weights back to G and every member of k_max at
809 # the maximum point. Also average the number of times that edge appears in
810 # the set of minimal 1-arborescences.
811 x_star = {}
812 size_k_max = len(k_max)
813 for u, v, d in G.edges(data=True):
814 edge_count = 0
815 d[weight] = original_edge_weights[(u, v)]
816 for k in k_max:
817 if (u, v) in k.edges():
818 edge_count += 1
819 k[u][v][weight] = original_edge_weights[(u, v)]
820 x_star[(u, v)] = edge_count / size_k_max
821 # Now symmetrize the edges in x_star and scale them according to (5) in
822 # reference [1]
823 z_star = {}
824 scale_factor = (G.order() - 1) / G.order()
825 for u, v in x_star:
826 frequency = x_star[(u, v)] + x_star[(v, u)]
827 if frequency > 0:
828 z_star[(u, v)] = scale_factor * frequency
829 del x_star
830 # Return the optimal weight and the z dict
831 return next(k_max.__iter__()).size(weight), z_star
834@nx._dispatchable
835def spanning_tree_distribution(G, z):
836 """
837 Find the asadpour exponential distribution of spanning trees.
839 Solves the Maximum Entropy Convex Program in the Asadpour algorithm [1]_
840 using the approach in section 7 to build an exponential distribution of
841 undirected spanning trees.
843 This algorithm ensures that the probability of any edge in a spanning
844 tree is proportional to the sum of the probabilities of the tress
845 containing that edge over the sum of the probabilities of all spanning
846 trees of the graph.
848 Parameters
849 ----------
850 G : nx.MultiGraph
851 The undirected support graph for the Held Karp relaxation
853 z : dict
854 The output of `held_karp_ascent()`, a scaled version of the Held-Karp
855 solution.
857 Returns
858 -------
859 gamma : dict
860 The probability distribution which approximately preserves the marginal
861 probabilities of `z`.
862 """
863 from math import exp
864 from math import log as ln
866 def q(e):
867 """
868 The value of q(e) is described in the Asadpour paper is "the
869 probability that edge e will be included in a spanning tree T that is
870 chosen with probability proportional to exp(gamma(T))" which
871 basically means that it is the total probability of the edge appearing
872 across the whole distribution.
874 Parameters
875 ----------
876 e : tuple
877 The `(u, v)` tuple describing the edge we are interested in
879 Returns
880 -------
881 float
882 The probability that a spanning tree chosen according to the
883 current values of gamma will include edge `e`.
884 """
885 # Create the laplacian matrices
886 for u, v, d in G.edges(data=True):
887 d[lambda_key] = exp(gamma[(u, v)])
888 G_Kirchhoff = nx.number_of_spanning_trees(G, weight=lambda_key)
889 G_e = nx.contracted_edge(G, e, self_loops=False)
890 G_e_Kirchhoff = nx.number_of_spanning_trees(G_e, weight=lambda_key)
892 # Multiply by the weight of the contracted edge since it is not included
893 # in the total weight of the contracted graph.
894 return exp(gamma[(e[0], e[1])]) * G_e_Kirchhoff / G_Kirchhoff
896 # initialize gamma to the zero dict
897 gamma = {}
898 for u, v, _ in G.edges:
899 gamma[(u, v)] = 0
901 # set epsilon
902 EPSILON = 0.2
904 # pick an edge attribute name that is unlikely to be in the graph
905 lambda_key = "spanning_tree_distribution's secret attribute name for lambda"
907 while True:
908 # We need to know that know that no values of q_e are greater than
909 # (1 + epsilon) * z_e, however changing one gamma value can increase the
910 # value of a different q_e, so we have to complete the for loop without
911 # changing anything for the condition to be meet
912 in_range_count = 0
913 # Search for an edge with q_e > (1 + epsilon) * z_e
914 for u, v in gamma:
915 e = (u, v)
916 q_e = q(e)
917 z_e = z[e]
918 if q_e > (1 + EPSILON) * z_e:
919 delta = ln(
920 (q_e * (1 - (1 + EPSILON / 2) * z_e))
921 / ((1 - q_e) * (1 + EPSILON / 2) * z_e)
922 )
923 gamma[e] -= delta
924 # Check that delta had the desired effect
925 new_q_e = q(e)
926 desired_q_e = (1 + EPSILON / 2) * z_e
927 if round(new_q_e, 8) != round(desired_q_e, 8):
928 raise nx.NetworkXError(
929 f"Unable to modify probability for edge ({u}, {v})"
930 )
931 else:
932 in_range_count += 1
933 # Check if the for loop terminated without changing any gamma
934 if in_range_count == len(gamma):
935 break
937 # Remove the new edge attributes
938 for _, _, d in G.edges(data=True):
939 if lambda_key in d:
940 del d[lambda_key]
942 return gamma
945@nx._dispatchable(edge_attrs="weight")
946def greedy_tsp(G, weight="weight", source=None):
947 """Return a low cost cycle starting at `source` and its cost.
949 This approximates a solution to the traveling salesman problem.
950 It finds a cycle of all the nodes that a salesman can visit in order
951 to visit many nodes while minimizing total distance.
952 It uses a simple greedy algorithm.
953 In essence, this function returns a large cycle given a source point
954 for which the total cost of the cycle is minimized.
956 Parameters
957 ----------
958 G : Graph
959 The Graph should be a complete weighted undirected graph.
960 The distance between all pairs of nodes should be included.
962 weight : string, optional (default="weight")
963 Edge data key corresponding to the edge weight.
964 If any edge does not have this attribute the weight is set to 1.
966 source : node, optional (default: first node in list(G))
967 Starting node. If None, defaults to ``next(iter(G))``
969 Returns
970 -------
971 cycle : list of nodes
972 Returns the cycle (list of nodes) that a salesman
973 can follow to minimize total weight of the trip.
975 Raises
976 ------
977 NetworkXError
978 If `G` is not complete, the algorithm raises an exception.
980 Examples
981 --------
982 >>> from networkx.algorithms import approximation as approx
983 >>> G = nx.DiGraph()
984 >>> G.add_weighted_edges_from(
985 ... {
986 ... ("A", "B", 3),
987 ... ("A", "C", 17),
988 ... ("A", "D", 14),
989 ... ("B", "A", 3),
990 ... ("B", "C", 12),
991 ... ("B", "D", 16),
992 ... ("C", "A", 13),
993 ... ("C", "B", 12),
994 ... ("C", "D", 4),
995 ... ("D", "A", 14),
996 ... ("D", "B", 15),
997 ... ("D", "C", 2),
998 ... }
999 ... )
1000 >>> cycle = approx.greedy_tsp(G, source="D")
1001 >>> cost = sum(G[n][nbr]["weight"] for n, nbr in nx.utils.pairwise(cycle))
1002 >>> cycle
1003 ['D', 'C', 'B', 'A', 'D']
1004 >>> cost
1005 31
1007 Notes
1008 -----
1009 This implementation of a greedy algorithm is based on the following:
1011 - The algorithm adds a node to the solution at every iteration.
1012 - The algorithm selects a node not already in the cycle whose connection
1013 to the previous node adds the least cost to the cycle.
1015 A greedy algorithm does not always give the best solution.
1016 However, it can construct a first feasible solution which can
1017 be passed as a parameter to an iterative improvement algorithm such
1018 as Simulated Annealing, or Threshold Accepting.
1020 Time complexity: It has a running time $O(|V|^2)$
1021 """
1022 # Check that G is a complete graph
1023 N = len(G) - 1
1024 # This check ignores selfloops which is what we want here.
1025 if any(len(nbrdict) - (n in nbrdict) != N for n, nbrdict in G.adj.items()):
1026 raise nx.NetworkXError("G must be a complete graph.")
1028 if source is None:
1029 source = nx.utils.arbitrary_element(G)
1031 if G.number_of_nodes() == 2:
1032 neighbor = next(G.neighbors(source))
1033 return [source, neighbor, source]
1035 nodeset = set(G)
1036 nodeset.remove(source)
1037 cycle = [source]
1038 next_node = source
1039 while nodeset:
1040 nbrdict = G[next_node]
1041 next_node = min(nodeset, key=lambda n: nbrdict[n].get(weight, 1))
1042 cycle.append(next_node)
1043 nodeset.remove(next_node)
1044 cycle.append(cycle[0])
1045 return cycle
1048@py_random_state(9)
1049@nx._dispatchable(edge_attrs="weight")
1050def simulated_annealing_tsp(
1051 G,
1052 init_cycle,
1053 weight="weight",
1054 source=None,
1055 temp=100,
1056 move="1-1",
1057 max_iterations=10,
1058 N_inner=100,
1059 alpha=0.01,
1060 seed=None,
1061):
1062 """Returns an approximate solution to the traveling salesman problem.
1064 This function uses simulated annealing to approximate the minimal cost
1065 cycle through the nodes. Starting from a suboptimal solution, simulated
1066 annealing perturbs that solution, occasionally accepting changes that make
1067 the solution worse to escape from a locally optimal solution. The chance
1068 of accepting such changes decreases over the iterations to encourage
1069 an optimal result. In summary, the function returns a cycle starting
1070 at `source` for which the total cost is minimized. It also returns the cost.
1072 The chance of accepting a proposed change is related to a parameter called
1073 the temperature (annealing has a physical analogue of steel hardening
1074 as it cools). As the temperature is reduced, the chance of moves that
1075 increase cost goes down.
1077 Parameters
1078 ----------
1079 G : Graph
1080 `G` should be a complete weighted graph.
1081 The distance between all pairs of nodes should be included.
1083 init_cycle : list of all nodes or "greedy"
1084 The initial solution (a cycle through all nodes returning to the start).
1085 This argument has no default to make you think about it.
1086 If "greedy", use `greedy_tsp(G, weight)`.
1087 Other common starting cycles are `list(G) + [next(iter(G))]` or the final
1088 result of `simulated_annealing_tsp` when doing `threshold_accepting_tsp`.
1090 weight : string, optional (default="weight")
1091 Edge data key corresponding to the edge weight.
1092 If any edge does not have this attribute the weight is set to 1.
1094 source : node, optional (default: first node in list(G))
1095 Starting node. If None, defaults to ``next(iter(G))``
1097 temp : int, optional (default=100)
1098 The algorithm's temperature parameter. It represents the initial
1099 value of temperature
1101 move : "1-1" or "1-0" or function, optional (default="1-1")
1102 Indicator of what move to use when finding new trial solutions.
1103 Strings indicate two special built-in moves:
1105 - "1-1": 1-1 exchange which transposes the position
1106 of two elements of the current solution.
1107 The function called is :func:`swap_two_nodes`.
1108 For example if we apply 1-1 exchange in the solution
1109 ``A = [3, 2, 1, 4, 3]``
1110 we can get the following by the transposition of 1 and 4 elements:
1111 ``A' = [3, 2, 4, 1, 3]``
1112 - "1-0": 1-0 exchange which moves an node in the solution
1113 to a new position.
1114 The function called is :func:`move_one_node`.
1115 For example if we apply 1-0 exchange in the solution
1116 ``A = [3, 2, 1, 4, 3]``
1117 we can transfer the fourth element to the second position:
1118 ``A' = [3, 4, 2, 1, 3]``
1120 You may provide your own functions to enact a move from
1121 one solution to a neighbor solution. The function must take
1122 the solution as input along with a `seed` input to control
1123 random number generation (see the `seed` input here).
1124 Your function should maintain the solution as a cycle with
1125 equal first and last node and all others appearing once.
1126 Your function should return the new solution.
1128 max_iterations : int, optional (default=10)
1129 Declared done when this number of consecutive iterations of
1130 the outer loop occurs without any change in the best cost solution.
1132 N_inner : int, optional (default=100)
1133 The number of iterations of the inner loop.
1135 alpha : float between (0, 1), optional (default=0.01)
1136 Percentage of temperature decrease in each iteration
1137 of outer loop
1139 seed : integer, random_state, or None (default)
1140 Indicator of random number generation state.
1141 See :ref:`Randomness<randomness>`.
1143 Returns
1144 -------
1145 cycle : list of nodes
1146 Returns the cycle (list of nodes) that a salesman
1147 can follow to minimize total weight of the trip.
1149 Raises
1150 ------
1151 NetworkXError
1152 If `G` is not complete the algorithm raises an exception.
1154 Examples
1155 --------
1156 >>> from networkx.algorithms import approximation as approx
1157 >>> G = nx.DiGraph()
1158 >>> G.add_weighted_edges_from(
1159 ... {
1160 ... ("A", "B", 3),
1161 ... ("A", "C", 17),
1162 ... ("A", "D", 14),
1163 ... ("B", "A", 3),
1164 ... ("B", "C", 12),
1165 ... ("B", "D", 16),
1166 ... ("C", "A", 13),
1167 ... ("C", "B", 12),
1168 ... ("C", "D", 4),
1169 ... ("D", "A", 14),
1170 ... ("D", "B", 15),
1171 ... ("D", "C", 2),
1172 ... }
1173 ... )
1174 >>> cycle = approx.simulated_annealing_tsp(G, "greedy", source="D")
1175 >>> cost = sum(G[n][nbr]["weight"] for n, nbr in nx.utils.pairwise(cycle))
1176 >>> cycle
1177 ['D', 'C', 'B', 'A', 'D']
1178 >>> cost
1179 31
1180 >>> incycle = ["D", "B", "A", "C", "D"]
1181 >>> cycle = approx.simulated_annealing_tsp(G, incycle, source="D")
1182 >>> cost = sum(G[n][nbr]["weight"] for n, nbr in nx.utils.pairwise(cycle))
1183 >>> cycle
1184 ['D', 'C', 'B', 'A', 'D']
1185 >>> cost
1186 31
1188 Notes
1189 -----
1190 Simulated Annealing is a metaheuristic local search algorithm.
1191 The main characteristic of this algorithm is that it accepts
1192 even solutions which lead to the increase of the cost in order
1193 to escape from low quality local optimal solutions.
1195 This algorithm needs an initial solution. If not provided, it is
1196 constructed by a simple greedy algorithm. At every iteration, the
1197 algorithm selects thoughtfully a neighbor solution.
1198 Consider $c(x)$ cost of current solution and $c(x')$ cost of a
1199 neighbor solution.
1200 If $c(x') - c(x) <= 0$ then the neighbor solution becomes the current
1201 solution for the next iteration. Otherwise, the algorithm accepts
1202 the neighbor solution with probability $p = exp - ([c(x') - c(x)] / temp)$.
1203 Otherwise the current solution is retained.
1205 `temp` is a parameter of the algorithm and represents temperature.
1207 Time complexity:
1208 For $N_i$ iterations of the inner loop and $N_o$ iterations of the
1209 outer loop, this algorithm has running time $O(N_i * N_o * |V|)$.
1211 For more information and how the algorithm is inspired see:
1212 http://en.wikipedia.org/wiki/Simulated_annealing
1213 """
1214 if move == "1-1":
1215 move = swap_two_nodes
1216 elif move == "1-0":
1217 move = move_one_node
1218 if init_cycle == "greedy":
1219 # Construct an initial solution using a greedy algorithm.
1220 cycle = greedy_tsp(G, weight=weight, source=source)
1221 if G.number_of_nodes() == 2:
1222 return cycle
1224 else:
1225 cycle = list(init_cycle)
1226 if source is None:
1227 source = cycle[0]
1228 elif source != cycle[0]:
1229 raise nx.NetworkXError("source must be first node in init_cycle")
1230 if cycle[0] != cycle[-1]:
1231 raise nx.NetworkXError("init_cycle must be a cycle. (return to start)")
1233 if len(cycle) - 1 != len(G) or len(set(G.nbunch_iter(cycle))) != len(G):
1234 raise nx.NetworkXError("init_cycle should be a cycle over all nodes in G.")
1236 # Check that G is a complete graph
1237 N = len(G) - 1
1238 # This check ignores selfloops which is what we want here.
1239 if any(len(nbrdict) - (n in nbrdict) != N for n, nbrdict in G.adj.items()):
1240 raise nx.NetworkXError("G must be a complete graph.")
1242 if G.number_of_nodes() == 2:
1243 neighbor = next(G.neighbors(source))
1244 return [source, neighbor, source]
1246 # Find the cost of initial solution
1247 cost = sum(G[u][v].get(weight, 1) for u, v in pairwise(cycle))
1249 count = 0
1250 best_cycle = cycle.copy()
1251 best_cost = cost
1252 while count <= max_iterations and temp > 0:
1253 count += 1
1254 for i in range(N_inner):
1255 adj_sol = move(cycle, seed)
1256 adj_cost = sum(G[u][v].get(weight, 1) for u, v in pairwise(adj_sol))
1257 delta = adj_cost - cost
1258 if delta <= 0:
1259 # Set current solution the adjacent solution.
1260 cycle = adj_sol
1261 cost = adj_cost
1263 if cost < best_cost:
1264 count = 0
1265 best_cycle = cycle.copy()
1266 best_cost = cost
1267 else:
1268 # Accept even a worse solution with probability p.
1269 p = math.exp(-delta / temp)
1270 if p >= seed.random():
1271 cycle = adj_sol
1272 cost = adj_cost
1273 temp -= temp * alpha
1275 return best_cycle
1278@py_random_state(9)
1279@nx._dispatchable(edge_attrs="weight")
1280def threshold_accepting_tsp(
1281 G,
1282 init_cycle,
1283 weight="weight",
1284 source=None,
1285 threshold=1,
1286 move="1-1",
1287 max_iterations=10,
1288 N_inner=100,
1289 alpha=0.1,
1290 seed=None,
1291):
1292 """Returns an approximate solution to the traveling salesman problem.
1294 This function uses threshold accepting methods to approximate the minimal cost
1295 cycle through the nodes. Starting from a suboptimal solution, threshold
1296 accepting methods perturb that solution, accepting any changes that make
1297 the solution no worse than increasing by a threshold amount. Improvements
1298 in cost are accepted, but so are changes leading to small increases in cost.
1299 This allows the solution to leave suboptimal local minima in solution space.
1300 The threshold is decreased slowly as iterations proceed helping to ensure
1301 an optimum. In summary, the function returns a cycle starting at `source`
1302 for which the total cost is minimized.
1304 Parameters
1305 ----------
1306 G : Graph
1307 `G` should be a complete weighted graph.
1308 The distance between all pairs of nodes should be included.
1310 init_cycle : list or "greedy"
1311 The initial solution (a cycle through all nodes returning to the start).
1312 This argument has no default to make you think about it.
1313 If "greedy", use `greedy_tsp(G, weight)`.
1314 Other common starting cycles are `list(G) + [next(iter(G))]` or the final
1315 result of `simulated_annealing_tsp` when doing `threshold_accepting_tsp`.
1317 weight : string, optional (default="weight")
1318 Edge data key corresponding to the edge weight.
1319 If any edge does not have this attribute the weight is set to 1.
1321 source : node, optional (default: first node in list(G))
1322 Starting node. If None, defaults to ``next(iter(G))``
1324 threshold : int, optional (default=1)
1325 The algorithm's threshold parameter. It represents the initial
1326 threshold's value
1328 move : "1-1" or "1-0" or function, optional (default="1-1")
1329 Indicator of what move to use when finding new trial solutions.
1330 Strings indicate two special built-in moves:
1332 - "1-1": 1-1 exchange which transposes the position
1333 of two elements of the current solution.
1334 The function called is :func:`swap_two_nodes`.
1335 For example if we apply 1-1 exchange in the solution
1336 ``A = [3, 2, 1, 4, 3]``
1337 we can get the following by the transposition of 1 and 4 elements:
1338 ``A' = [3, 2, 4, 1, 3]``
1339 - "1-0": 1-0 exchange which moves an node in the solution
1340 to a new position.
1341 The function called is :func:`move_one_node`.
1342 For example if we apply 1-0 exchange in the solution
1343 ``A = [3, 2, 1, 4, 3]``
1344 we can transfer the fourth element to the second position:
1345 ``A' = [3, 4, 2, 1, 3]``
1347 You may provide your own functions to enact a move from
1348 one solution to a neighbor solution. The function must take
1349 the solution as input along with a `seed` input to control
1350 random number generation (see the `seed` input here).
1351 Your function should maintain the solution as a cycle with
1352 equal first and last node and all others appearing once.
1353 Your function should return the new solution.
1355 max_iterations : int, optional (default=10)
1356 Declared done when this number of consecutive iterations of
1357 the outer loop occurs without any change in the best cost solution.
1359 N_inner : int, optional (default=100)
1360 The number of iterations of the inner loop.
1362 alpha : float between (0, 1), optional (default=0.1)
1363 Percentage of threshold decrease when there is at
1364 least one acceptance of a neighbor solution.
1365 If no inner loop moves are accepted the threshold remains unchanged.
1367 seed : integer, random_state, or None (default)
1368 Indicator of random number generation state.
1369 See :ref:`Randomness<randomness>`.
1371 Returns
1372 -------
1373 cycle : list of nodes
1374 Returns the cycle (list of nodes) that a salesman
1375 can follow to minimize total weight of the trip.
1377 Raises
1378 ------
1379 NetworkXError
1380 If `G` is not complete the algorithm raises an exception.
1382 Examples
1383 --------
1384 >>> from networkx.algorithms import approximation as approx
1385 >>> G = nx.DiGraph()
1386 >>> G.add_weighted_edges_from(
1387 ... {
1388 ... ("A", "B", 3),
1389 ... ("A", "C", 17),
1390 ... ("A", "D", 14),
1391 ... ("B", "A", 3),
1392 ... ("B", "C", 12),
1393 ... ("B", "D", 16),
1394 ... ("C", "A", 13),
1395 ... ("C", "B", 12),
1396 ... ("C", "D", 4),
1397 ... ("D", "A", 14),
1398 ... ("D", "B", 15),
1399 ... ("D", "C", 2),
1400 ... }
1401 ... )
1402 >>> cycle = approx.threshold_accepting_tsp(G, "greedy", source="D")
1403 >>> cost = sum(G[n][nbr]["weight"] for n, nbr in nx.utils.pairwise(cycle))
1404 >>> cycle
1405 ['D', 'C', 'B', 'A', 'D']
1406 >>> cost
1407 31
1408 >>> incycle = ["D", "B", "A", "C", "D"]
1409 >>> cycle = approx.threshold_accepting_tsp(G, incycle, source="D")
1410 >>> cost = sum(G[n][nbr]["weight"] for n, nbr in nx.utils.pairwise(cycle))
1411 >>> cycle
1412 ['D', 'C', 'B', 'A', 'D']
1413 >>> cost
1414 31
1416 Notes
1417 -----
1418 Threshold Accepting is a metaheuristic local search algorithm.
1419 The main characteristic of this algorithm is that it accepts
1420 even solutions which lead to the increase of the cost in order
1421 to escape from low quality local optimal solutions.
1423 This algorithm needs an initial solution. This solution can be
1424 constructed by a simple greedy algorithm. At every iteration, it
1425 selects thoughtfully a neighbor solution.
1426 Consider $c(x)$ cost of current solution and $c(x')$ cost of
1427 neighbor solution.
1428 If $c(x') - c(x) <= threshold$ then the neighbor solution becomes the current
1429 solution for the next iteration, where the threshold is named threshold.
1431 In comparison to the Simulated Annealing algorithm, the Threshold
1432 Accepting algorithm does not accept very low quality solutions
1433 (due to the presence of the threshold value). In the case of
1434 Simulated Annealing, even a very low quality solution can
1435 be accepted with probability $p$.
1437 Time complexity:
1438 It has a running time $O(m * n * |V|)$ where $m$ and $n$ are the number
1439 of times the outer and inner loop run respectively.
1441 For more information and how algorithm is inspired see:
1442 https://doi.org/10.1016/0021-9991(90)90201-B
1444 See Also
1445 --------
1446 simulated_annealing_tsp
1448 """
1449 if move == "1-1":
1450 move = swap_two_nodes
1451 elif move == "1-0":
1452 move = move_one_node
1453 if init_cycle == "greedy":
1454 # Construct an initial solution using a greedy algorithm.
1455 cycle = greedy_tsp(G, weight=weight, source=source)
1456 if G.number_of_nodes() == 2:
1457 return cycle
1459 else:
1460 cycle = list(init_cycle)
1461 if source is None:
1462 source = cycle[0]
1463 elif source != cycle[0]:
1464 raise nx.NetworkXError("source must be first node in init_cycle")
1465 if cycle[0] != cycle[-1]:
1466 raise nx.NetworkXError("init_cycle must be a cycle. (return to start)")
1468 if len(cycle) - 1 != len(G) or len(set(G.nbunch_iter(cycle))) != len(G):
1469 raise nx.NetworkXError("init_cycle is not all and only nodes.")
1471 # Check that G is a complete graph
1472 N = len(G) - 1
1473 # This check ignores selfloops which is what we want here.
1474 if any(len(nbrdict) - (n in nbrdict) != N for n, nbrdict in G.adj.items()):
1475 raise nx.NetworkXError("G must be a complete graph.")
1477 if G.number_of_nodes() == 2:
1478 neighbor = list(G.neighbors(source))[0]
1479 return [source, neighbor, source]
1481 # Find the cost of initial solution
1482 cost = sum(G[u][v].get(weight, 1) for u, v in pairwise(cycle))
1484 count = 0
1485 best_cycle = cycle.copy()
1486 best_cost = cost
1487 while count <= max_iterations:
1488 count += 1
1489 accepted = False
1490 for i in range(N_inner):
1491 adj_sol = move(cycle, seed)
1492 adj_cost = sum(G[u][v].get(weight, 1) for u, v in pairwise(adj_sol))
1493 delta = adj_cost - cost
1494 if delta <= threshold:
1495 accepted = True
1497 # Set current solution the adjacent solution.
1498 cycle = adj_sol
1499 cost = adj_cost
1501 if cost < best_cost:
1502 count = 0
1503 best_cycle = cycle.copy()
1504 best_cost = cost
1505 if accepted:
1506 threshold -= threshold * alpha
1508 return best_cycle