Coverage for /pythoncovmergedfiles/medio/medio/usr/local/lib/python3.9/dist-packages/networkx/algorithms/approximation/traveling_salesman.py: 7%
420 statements
« prev ^ index » next coverage.py v7.3.2, created at 2023-10-20 07:00 +0000
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"""
36import math
38import networkx as nx
39from networkx.algorithms.tree.mst import random_spanning_tree
40from networkx.utils import not_implemented_for, pairwise, py_random_state
42__all__ = [
43 "traveling_salesman_problem",
44 "christofides",
45 "asadpour_atsp",
46 "greedy_tsp",
47 "simulated_annealing_tsp",
48 "threshold_accepting_tsp",
49]
52def swap_two_nodes(soln, seed):
53 """Swap two nodes in `soln` to give a neighbor solution.
55 Parameters
56 ----------
57 soln : list of nodes
58 Current cycle of nodes
60 seed : integer, random_state, or None (default)
61 Indicator of random number generation state.
62 See :ref:`Randomness<randomness>`.
64 Returns
65 -------
66 list
67 The solution after move is applied. (A neighbor solution.)
69 Notes
70 -----
71 This function assumes that the incoming list `soln` is a cycle
72 (that the first and last element are the same) and also that
73 we don't want any move to change the first node in the list
74 (and thus not the last node either).
76 The input list is changed as well as returned. Make a copy if needed.
78 See Also
79 --------
80 move_one_node
81 """
82 a, b = seed.sample(range(1, len(soln) - 1), k=2)
83 soln[a], soln[b] = soln[b], soln[a]
84 return soln
87def move_one_node(soln, seed):
88 """Move one node to another position to give a neighbor solution.
90 The node to move and the position to move to are chosen randomly.
91 The first and last nodes are left untouched as soln must be a cycle
92 starting at that node.
94 Parameters
95 ----------
96 soln : list of nodes
97 Current cycle of nodes
99 seed : integer, random_state, or None (default)
100 Indicator of random number generation state.
101 See :ref:`Randomness<randomness>`.
103 Returns
104 -------
105 list
106 The solution after move is applied. (A neighbor solution.)
108 Notes
109 -----
110 This function assumes that the incoming list `soln` is a cycle
111 (that the first and last element are the same) and also that
112 we don't want any move to change the first node in the list
113 (and thus not the last node either).
115 The input list is changed as well as returned. Make a copy if needed.
117 See Also
118 --------
119 swap_two_nodes
120 """
121 a, b = seed.sample(range(1, len(soln) - 1), k=2)
122 soln.insert(b, soln.pop(a))
123 return soln
126@not_implemented_for("directed")
127@nx._dispatch(edge_attrs="weight")
128def christofides(G, weight="weight", tree=None):
129 """Approximate a solution of the traveling salesman problem
131 Compute a 3/2-approximation of the traveling salesman problem
132 in a complete undirected graph using Christofides [1]_ algorithm.
134 Parameters
135 ----------
136 G : Graph
137 `G` should be a complete weighted undirected graph.
138 The distance between all pairs of nodes should be included.
140 weight : string, optional (default="weight")
141 Edge data key corresponding to the edge weight.
142 If any edge does not have this attribute the weight is set to 1.
144 tree : NetworkX graph or None (default: None)
145 A minimum spanning tree of G. Or, if None, the minimum spanning
146 tree is computed using :func:`networkx.minimum_spanning_tree`
148 Returns
149 -------
150 list
151 List of nodes in `G` along a cycle with a 3/2-approximation of
152 the minimal Hamiltonian cycle.
154 References
155 ----------
156 .. [1] Christofides, Nicos. "Worst-case analysis of a new heuristic for
157 the travelling salesman problem." No. RR-388. Carnegie-Mellon Univ
158 Pittsburgh Pa Management Sciences Research Group, 1976.
159 """
160 # Remove selfloops if necessary
161 loop_nodes = nx.nodes_with_selfloops(G)
162 try:
163 node = next(loop_nodes)
164 except StopIteration:
165 pass
166 else:
167 G = G.copy()
168 G.remove_edge(node, node)
169 G.remove_edges_from((n, n) for n in loop_nodes)
170 # Check that G is a complete graph
171 N = len(G) - 1
172 # This check ignores selfloops which is what we want here.
173 if any(len(nbrdict) != N for n, nbrdict in G.adj.items()):
174 raise nx.NetworkXError("G must be a complete graph.")
176 if tree is None:
177 tree = nx.minimum_spanning_tree(G, weight=weight)
178 L = G.copy()
179 L.remove_nodes_from([v for v, degree in tree.degree if not (degree % 2)])
180 MG = nx.MultiGraph()
181 MG.add_edges_from(tree.edges)
182 edges = nx.min_weight_matching(L, weight=weight)
183 MG.add_edges_from(edges)
184 return _shortcutting(nx.eulerian_circuit(MG))
187def _shortcutting(circuit):
188 """Remove duplicate nodes in the path"""
189 nodes = []
190 for u, v in circuit:
191 if v in nodes:
192 continue
193 if not nodes:
194 nodes.append(u)
195 nodes.append(v)
196 nodes.append(nodes[0])
197 return nodes
200@nx._dispatch(edge_attrs="weight")
201def traveling_salesman_problem(G, weight="weight", nodes=None, cycle=True, method=None):
202 """Find the shortest path in `G` connecting specified nodes
204 This function allows approximate solution to the traveling salesman
205 problem on networks that are not complete graphs and/or where the
206 salesman does not need to visit all nodes.
208 This function proceeds in two steps. First, it creates a complete
209 graph using the all-pairs shortest_paths between nodes in `nodes`.
210 Edge weights in the new graph are the lengths of the paths
211 between each pair of nodes in the original graph.
212 Second, an algorithm (default: `christofides` for undirected and
213 `asadpour_atsp` for directed) is used to approximate the minimal Hamiltonian
214 cycle on this new graph. The available algorithms are:
216 - christofides
217 - greedy_tsp
218 - simulated_annealing_tsp
219 - threshold_accepting_tsp
220 - asadpour_atsp
222 Once the Hamiltonian Cycle is found, this function post-processes to
223 accommodate the structure of the original graph. If `cycle` is ``False``,
224 the biggest weight edge is removed to make a Hamiltonian path.
225 Then each edge on the new complete graph used for that analysis is
226 replaced by the shortest_path between those nodes on the original graph.
228 Parameters
229 ----------
230 G : NetworkX graph
231 A possibly weighted graph
233 nodes : collection of nodes (default=G.nodes)
234 collection (list, set, etc.) of nodes to visit
236 weight : string, optional (default="weight")
237 Edge data key corresponding to the edge weight.
238 If any edge does not have this attribute the weight is set to 1.
240 cycle : bool (default: True)
241 Indicates whether a cycle should be returned, or a path.
242 Note: the cycle is the approximate minimal cycle.
243 The path simply removes the biggest edge in that cycle.
245 method : function (default: None)
246 A function that returns a cycle on all nodes and approximates
247 the solution to the traveling salesman problem on a complete
248 graph. The returned cycle is then used to find a corresponding
249 solution on `G`. `method` should be callable; take inputs
250 `G`, and `weight`; and return a list of nodes along the cycle.
252 Provided options include :func:`christofides`, :func:`greedy_tsp`,
253 :func:`simulated_annealing_tsp` and :func:`threshold_accepting_tsp`.
255 If `method is None`: use :func:`christofides` for undirected `G` and
256 :func:`threshold_accepting_tsp` for directed `G`.
258 To specify parameters for these provided functions, construct lambda
259 functions that state the specific value. `method` must have 2 inputs.
260 (See examples).
262 Returns
263 -------
264 list
265 List of nodes in `G` along a path with an approximation of the minimal
266 path through `nodes`.
269 Raises
270 ------
271 NetworkXError
272 If `G` is a directed graph it has to be strongly connected or the
273 complete version cannot be generated.
275 Examples
276 --------
277 >>> tsp = nx.approximation.traveling_salesman_problem
278 >>> G = nx.cycle_graph(9)
279 >>> G[4][5]["weight"] = 5 # all other weights are 1
280 >>> tsp(G, nodes=[3, 6])
281 [3, 2, 1, 0, 8, 7, 6, 7, 8, 0, 1, 2, 3]
282 >>> path = tsp(G, cycle=False)
283 >>> path in ([4, 3, 2, 1, 0, 8, 7, 6, 5], [5, 6, 7, 8, 0, 1, 2, 3, 4])
284 True
286 Build (curry) your own function to provide parameter values to the methods.
288 >>> SA_tsp = nx.approximation.simulated_annealing_tsp
289 >>> method = lambda G, wt: SA_tsp(G, "greedy", weight=wt, temp=500)
290 >>> path = tsp(G, cycle=False, method=method)
291 >>> path in ([4, 3, 2, 1, 0, 8, 7, 6, 5], [5, 6, 7, 8, 0, 1, 2, 3, 4])
292 True
294 """
295 if method is None:
296 if G.is_directed():
297 method = asadpour_atsp
298 else:
299 method = christofides
300 if nodes is None:
301 nodes = list(G.nodes)
303 dist = {}
304 path = {}
305 for n, (d, p) in nx.all_pairs_dijkstra(G, weight=weight):
306 dist[n] = d
307 path[n] = p
309 if G.is_directed():
310 # If the graph is not strongly connected, raise an exception
311 if not nx.is_strongly_connected(G):
312 raise nx.NetworkXError("G is not strongly connected")
313 GG = nx.DiGraph()
314 else:
315 GG = nx.Graph()
316 for u in nodes:
317 for v in nodes:
318 if u == v:
319 continue
320 GG.add_edge(u, v, weight=dist[u][v])
321 best_GG = method(GG, weight)
323 if not cycle:
324 # find and remove the biggest edge
325 (u, v) = max(pairwise(best_GG), key=lambda x: dist[x[0]][x[1]])
326 pos = best_GG.index(u) + 1
327 while best_GG[pos] != v:
328 pos = best_GG[pos:].index(u) + 1
329 best_GG = best_GG[pos:-1] + best_GG[:pos]
331 best_path = []
332 for u, v in pairwise(best_GG):
333 best_path.extend(path[u][v][:-1])
334 best_path.append(v)
335 return best_path
338@not_implemented_for("undirected")
339@py_random_state(2)
340@nx._dispatch(edge_attrs="weight")
341def asadpour_atsp(G, weight="weight", seed=None, source=None):
342 """
343 Returns an approximate solution to the traveling salesman problem.
345 This approximate solution is one of the best known approximations for the
346 asymmetric traveling salesman problem developed by Asadpour et al,
347 [1]_. The algorithm first solves the Held-Karp relaxation to find a lower
348 bound for the weight of the cycle. Next, it constructs an exponential
349 distribution of undirected spanning trees where the probability of an
350 edge being in the tree corresponds to the weight of that edge using a
351 maximum entropy rounding scheme. Next we sample that distribution
352 $2 \\lceil \\ln n \\rceil$ times and save the minimum sampled tree once the
353 direction of the arcs is added back to the edges. Finally, we augment
354 then short circuit that graph to find the approximate tour for the
355 salesman.
357 Parameters
358 ----------
359 G : nx.DiGraph
360 The graph should be a complete weighted directed graph. The
361 distance between all paris of nodes should be included and the triangle
362 inequality should hold. That is, the direct edge between any two nodes
363 should be the path of least cost.
365 weight : string, optional (default="weight")
366 Edge data key corresponding to the edge weight.
367 If any edge does not have this attribute the weight is set to 1.
369 seed : integer, random_state, or None (default)
370 Indicator of random number generation state.
371 See :ref:`Randomness<randomness>`.
373 source : node label (default=`None`)
374 If given, return the cycle starting and ending at the given node.
376 Returns
377 -------
378 cycle : list of nodes
379 Returns the cycle (list of nodes) that a salesman can follow to minimize
380 the total weight of the trip.
382 Raises
383 ------
384 NetworkXError
385 If `G` is not complete or has less than two nodes, the algorithm raises
386 an exception.
388 NetworkXError
389 If `source` is not `None` and is not a node in `G`, the algorithm raises
390 an exception.
392 NetworkXNotImplemented
393 If `G` is an undirected graph.
395 References
396 ----------
397 .. [1] A. Asadpour, M. X. Goemans, A. Madry, S. O. Gharan, and A. Saberi,
398 An o(log n/log log n)-approximation algorithm for the asymmetric
399 traveling salesman problem, Operations research, 65 (2017),
400 pp. 1043–1061
402 Examples
403 --------
404 >>> import networkx as nx
405 >>> import networkx.algorithms.approximation as approx
406 >>> G = nx.complete_graph(3, create_using=nx.DiGraph)
407 >>> nx.set_edge_attributes(G, {(0, 1): 2, (1, 2): 2, (2, 0): 2, (0, 2): 1, (2, 1): 1, (1, 0): 1}, "weight")
408 >>> tour = approx.asadpour_atsp(G,source=0)
409 >>> tour
410 [0, 2, 1, 0]
411 """
412 from math import ceil, exp
413 from math import log as ln
415 # Check that G is a complete graph
416 N = len(G) - 1
417 if N < 2:
418 raise nx.NetworkXError("G must have at least two nodes")
419 # This check ignores selfloops which is what we want here.
420 if any(len(nbrdict) - (n in nbrdict) != N for n, nbrdict in G.adj.items()):
421 raise nx.NetworkXError("G is not a complete DiGraph")
422 # Check that the source vertex, if given, is in the graph
423 if source is not None and source not in G.nodes:
424 raise nx.NetworkXError("Given source node not in G.")
426 opt_hk, z_star = held_karp_ascent(G, weight)
428 # Test to see if the ascent method found an integer solution or a fractional
429 # solution. If it is integral then z_star is a nx.Graph, otherwise it is
430 # a dict
431 if not isinstance(z_star, dict):
432 # Here we are using the shortcutting method to go from the list of edges
433 # returned from eulerian_circuit to a list of nodes
434 return _shortcutting(nx.eulerian_circuit(z_star, source=source))
436 # Create the undirected support of z_star
437 z_support = nx.MultiGraph()
438 for u, v in z_star:
439 if (u, v) not in z_support.edges:
440 edge_weight = min(G[u][v][weight], G[v][u][weight])
441 z_support.add_edge(u, v, **{weight: edge_weight})
443 # Create the exponential distribution of spanning trees
444 gamma = spanning_tree_distribution(z_support, z_star)
446 # Write the lambda values to the edges of z_support
447 z_support = nx.Graph(z_support)
448 lambda_dict = {(u, v): exp(gamma[(u, v)]) for u, v in z_support.edges()}
449 nx.set_edge_attributes(z_support, lambda_dict, "weight")
450 del gamma, lambda_dict
452 # Sample 2 * ceil( ln(n) ) spanning trees and record the minimum one
453 minimum_sampled_tree = None
454 minimum_sampled_tree_weight = math.inf
455 for _ in range(2 * ceil(ln(G.number_of_nodes()))):
456 sampled_tree = random_spanning_tree(z_support, "weight", seed=seed)
457 sampled_tree_weight = sampled_tree.size(weight)
458 if sampled_tree_weight < minimum_sampled_tree_weight:
459 minimum_sampled_tree = sampled_tree.copy()
460 minimum_sampled_tree_weight = sampled_tree_weight
462 # Orient the edges in that tree to keep the cost of the tree the same.
463 t_star = nx.MultiDiGraph()
464 for u, v, d in minimum_sampled_tree.edges(data=weight):
465 if d == G[u][v][weight]:
466 t_star.add_edge(u, v, **{weight: d})
467 else:
468 t_star.add_edge(v, u, **{weight: d})
470 # Find the node demands needed to neutralize the flow of t_star in G
471 node_demands = {n: t_star.out_degree(n) - t_star.in_degree(n) for n in t_star}
472 nx.set_node_attributes(G, node_demands, "demand")
474 # Find the min_cost_flow
475 flow_dict = nx.min_cost_flow(G, "demand")
477 # Build the flow into t_star
478 for source, values in flow_dict.items():
479 for target in values:
480 if (source, target) not in t_star.edges and values[target] > 0:
481 # IF values[target] > 0 we have to add that many edges
482 for _ in range(values[target]):
483 t_star.add_edge(source, target)
485 # Return the shortcut eulerian circuit
486 circuit = nx.eulerian_circuit(t_star, source=source)
487 return _shortcutting(circuit)
490@nx._dispatch(edge_attrs="weight")
491def held_karp_ascent(G, weight="weight"):
492 """
493 Minimizes the Held-Karp relaxation of the TSP for `G`
495 Solves the Held-Karp relaxation of the input complete digraph and scales
496 the output solution for use in the Asadpour [1]_ ASTP algorithm.
498 The Held-Karp relaxation defines the lower bound for solutions to the
499 ATSP, although it does return a fractional solution. This is used in the
500 Asadpour algorithm as an initial solution which is later rounded to a
501 integral tree within the spanning tree polytopes. This function solves
502 the relaxation with the branch and bound method in [2]_.
504 Parameters
505 ----------
506 G : nx.DiGraph
507 The graph should be a complete weighted directed graph.
508 The distance between all paris of nodes should be included.
510 weight : string, optional (default="weight")
511 Edge data key corresponding to the edge weight.
512 If any edge does not have this attribute the weight is set to 1.
514 Returns
515 -------
516 OPT : float
517 The cost for the optimal solution to the Held-Karp relaxation
518 z : dict or nx.Graph
519 A symmetrized and scaled version of the optimal solution to the
520 Held-Karp relaxation for use in the Asadpour algorithm.
522 If an integral solution is found, then that is an optimal solution for
523 the ATSP problem and that is returned instead.
525 References
526 ----------
527 .. [1] A. Asadpour, M. X. Goemans, A. Madry, S. O. Gharan, and A. Saberi,
528 An o(log n/log log n)-approximation algorithm for the asymmetric
529 traveling salesman problem, Operations research, 65 (2017),
530 pp. 1043–1061
532 .. [2] M. Held, R. M. Karp, The traveling-salesman problem and minimum
533 spanning trees, Operations Research, 1970-11-01, Vol. 18 (6),
534 pp.1138-1162
535 """
536 import numpy as np
537 from scipy import optimize
539 def k_pi():
540 """
541 Find the set of minimum 1-Arborescences for G at point pi.
543 Returns
544 -------
545 Set
546 The set of minimum 1-Arborescences
547 """
548 # Create a copy of G without vertex 1.
549 G_1 = G.copy()
550 minimum_1_arborescences = set()
551 minimum_1_arborescence_weight = math.inf
553 # node is node '1' in the Held and Karp paper
554 n = next(G.__iter__())
555 G_1.remove_node(n)
557 # Iterate over the spanning arborescences of the graph until we know
558 # that we have found the minimum 1-arborescences. My proposed strategy
559 # is to find the most extensive root to connect to from 'node 1' and
560 # the least expensive one. We then iterate over arborescences until
561 # the cost of the basic arborescence is the cost of the minimum one
562 # plus the difference between the most and least expensive roots,
563 # that way the cost of connecting 'node 1' will by definition not by
564 # minimum
565 min_root = {"node": None, weight: math.inf}
566 max_root = {"node": None, weight: -math.inf}
567 for u, v, d in G.edges(n, data=True):
568 if d[weight] < min_root[weight]:
569 min_root = {"node": v, weight: d[weight]}
570 if d[weight] > max_root[weight]:
571 max_root = {"node": v, weight: d[weight]}
573 min_in_edge = min(G.in_edges(n, data=True), key=lambda x: x[2][weight])
574 min_root[weight] = min_root[weight] + min_in_edge[2][weight]
575 max_root[weight] = max_root[weight] + min_in_edge[2][weight]
577 min_arb_weight = math.inf
578 for arb in nx.ArborescenceIterator(G_1):
579 arb_weight = arb.size(weight)
580 if min_arb_weight == math.inf:
581 min_arb_weight = arb_weight
582 elif arb_weight > min_arb_weight + max_root[weight] - min_root[weight]:
583 break
584 # We have to pick the root node of the arborescence for the out
585 # edge of the first vertex as that is the only node without an
586 # edge directed into it.
587 for N, deg in arb.in_degree:
588 if deg == 0:
589 # root found
590 arb.add_edge(n, N, **{weight: G[n][N][weight]})
591 arb_weight += G[n][N][weight]
592 break
594 # We can pick the minimum weight in-edge for the vertex with
595 # a cycle. If there are multiple edges with the same, minimum
596 # weight, We need to add all of them.
597 #
598 # Delete the edge (N, v) so that we cannot pick it.
599 edge_data = G[N][n]
600 G.remove_edge(N, n)
601 min_weight = min(G.in_edges(n, data=weight), key=lambda x: x[2])[2]
602 min_edges = [
603 (u, v, d) for u, v, d in G.in_edges(n, data=weight) if d == min_weight
604 ]
605 for u, v, d in min_edges:
606 new_arb = arb.copy()
607 new_arb.add_edge(u, v, **{weight: d})
608 new_arb_weight = arb_weight + d
609 # Check to see the weight of the arborescence, if it is a
610 # new minimum, clear all of the old potential minimum
611 # 1-arborescences and add this is the only one. If its
612 # weight is above the known minimum, do not add it.
613 if new_arb_weight < minimum_1_arborescence_weight:
614 minimum_1_arborescences.clear()
615 minimum_1_arborescence_weight = new_arb_weight
616 # We have a 1-arborescence, add it to the set
617 if new_arb_weight == minimum_1_arborescence_weight:
618 minimum_1_arborescences.add(new_arb)
619 G.add_edge(N, n, **edge_data)
621 return minimum_1_arborescences
623 def direction_of_ascent():
624 """
625 Find the direction of ascent at point pi.
627 See [1]_ for more information.
629 Returns
630 -------
631 dict
632 A mapping from the nodes of the graph which represents the direction
633 of ascent.
635 References
636 ----------
637 .. [1] M. Held, R. M. Karp, The traveling-salesman problem and minimum
638 spanning trees, Operations Research, 1970-11-01, Vol. 18 (6),
639 pp.1138-1162
640 """
641 # 1. Set d equal to the zero n-vector.
642 d = {}
643 for n in G:
644 d[n] = 0
645 del n
646 # 2. Find a 1-Arborescence T^k such that k is in K(pi, d).
647 minimum_1_arborescences = k_pi()
648 while True:
649 # Reduce K(pi) to K(pi, d)
650 # Find the arborescence in K(pi) which increases the lest in
651 # direction d
652 min_k_d_weight = math.inf
653 min_k_d = None
654 for arborescence in minimum_1_arborescences:
655 weighted_cost = 0
656 for n, deg in arborescence.degree:
657 weighted_cost += d[n] * (deg - 2)
658 if weighted_cost < min_k_d_weight:
659 min_k_d_weight = weighted_cost
660 min_k_d = arborescence
662 # 3. If sum of d_i * v_{i, k} is greater than zero, terminate
663 if min_k_d_weight > 0:
664 return d, min_k_d
665 # 4. d_i = d_i + v_{i, k}
666 for n, deg in min_k_d.degree:
667 d[n] += deg - 2
668 # Check that we do not need to terminate because the direction
669 # of ascent does not exist. This is done with linear
670 # programming.
671 c = np.full(len(minimum_1_arborescences), -1, dtype=int)
672 a_eq = np.empty((len(G) + 1, len(minimum_1_arborescences)), dtype=int)
673 b_eq = np.zeros(len(G) + 1, dtype=int)
674 b_eq[len(G)] = 1
675 for arb_count, arborescence in enumerate(minimum_1_arborescences):
676 n_count = len(G) - 1
677 for n, deg in arborescence.degree:
678 a_eq[n_count][arb_count] = deg - 2
679 n_count -= 1
680 a_eq[len(G)][arb_count] = 1
681 program_result = optimize.linprog(c, A_eq=a_eq, b_eq=b_eq)
682 # If the constants exist, then the direction of ascent doesn't
683 if program_result.success:
684 # There is no direction of ascent
685 return None, minimum_1_arborescences
687 # 5. GO TO 2
689 def find_epsilon(k, d):
690 """
691 Given the direction of ascent at pi, find the maximum distance we can go
692 in that direction.
694 Parameters
695 ----------
696 k_xy : set
697 The set of 1-arborescences which have the minimum rate of increase
698 in the direction of ascent
700 d : dict
701 The direction of ascent
703 Returns
704 -------
705 float
706 The distance we can travel in direction `d`
707 """
708 min_epsilon = math.inf
709 for e_u, e_v, e_w in G.edges(data=weight):
710 if (e_u, e_v) in k.edges:
711 continue
712 # Now, I have found a condition which MUST be true for the edges to
713 # be a valid substitute. The edge in the graph which is the
714 # substitute is the one with the same terminated end. This can be
715 # checked rather simply.
716 #
717 # Find the edge within k which is the substitute. Because k is a
718 # 1-arborescence, we know that they is only one such edges
719 # leading into every vertex.
720 if len(k.in_edges(e_v, data=weight)) > 1:
721 raise Exception
722 sub_u, sub_v, sub_w = next(k.in_edges(e_v, data=weight).__iter__())
723 k.add_edge(e_u, e_v, **{weight: e_w})
724 k.remove_edge(sub_u, sub_v)
725 if (
726 max(d for n, d in k.in_degree()) <= 1
727 and len(G) == k.number_of_edges()
728 and nx.is_weakly_connected(k)
729 ):
730 # Ascent method calculation
731 if d[sub_u] == d[e_u] or sub_w == e_w:
732 # Revert to the original graph
733 k.remove_edge(e_u, e_v)
734 k.add_edge(sub_u, sub_v, **{weight: sub_w})
735 continue
736 epsilon = (sub_w - e_w) / (d[e_u] - d[sub_u])
737 if 0 < epsilon < min_epsilon:
738 min_epsilon = epsilon
739 # Revert to the original graph
740 k.remove_edge(e_u, e_v)
741 k.add_edge(sub_u, sub_v, **{weight: sub_w})
743 return min_epsilon
745 # I have to know that the elements in pi correspond to the correct elements
746 # in the direction of ascent, even if the node labels are not integers.
747 # Thus, I will use dictionaries to made that mapping.
748 pi_dict = {}
749 for n in G:
750 pi_dict[n] = 0
751 del n
752 original_edge_weights = {}
753 for u, v, d in G.edges(data=True):
754 original_edge_weights[(u, v)] = d[weight]
755 dir_ascent, k_d = direction_of_ascent()
756 while dir_ascent is not None:
757 max_distance = find_epsilon(k_d, dir_ascent)
758 for n, v in dir_ascent.items():
759 pi_dict[n] += max_distance * v
760 for u, v, d in G.edges(data=True):
761 d[weight] = original_edge_weights[(u, v)] + pi_dict[u]
762 dir_ascent, k_d = direction_of_ascent()
763 # k_d is no longer an individual 1-arborescence but rather a set of
764 # minimal 1-arborescences at the maximum point of the polytope and should
765 # be reflected as such
766 k_max = k_d
768 # Search for a cycle within k_max. If a cycle exists, return it as the
769 # solution
770 for k in k_max:
771 if len([n for n in k if k.degree(n) == 2]) == G.order():
772 # Tour found
773 return k.size(weight), k
775 # Write the original edge weights back to G and every member of k_max at
776 # the maximum point. Also average the number of times that edge appears in
777 # the set of minimal 1-arborescences.
778 x_star = {}
779 size_k_max = len(k_max)
780 for u, v, d in G.edges(data=True):
781 edge_count = 0
782 d[weight] = original_edge_weights[(u, v)]
783 for k in k_max:
784 if (u, v) in k.edges():
785 edge_count += 1
786 k[u][v][weight] = original_edge_weights[(u, v)]
787 x_star[(u, v)] = edge_count / size_k_max
788 # Now symmetrize the edges in x_star and scale them according to (5) in
789 # reference [1]
790 z_star = {}
791 scale_factor = (G.order() - 1) / G.order()
792 for u, v in x_star:
793 frequency = x_star[(u, v)] + x_star[(v, u)]
794 if frequency > 0:
795 z_star[(u, v)] = scale_factor * frequency
796 del x_star
797 # Return the optimal weight and the z dict
798 return next(k_max.__iter__()).size(weight), z_star
801@nx._dispatch
802def spanning_tree_distribution(G, z):
803 """
804 Find the asadpour exponential distribution of spanning trees.
806 Solves the Maximum Entropy Convex Program in the Asadpour algorithm [1]_
807 using the approach in section 7 to build an exponential distribution of
808 undirected spanning trees.
810 This algorithm ensures that the probability of any edge in a spanning
811 tree is proportional to the sum of the probabilities of the tress
812 containing that edge over the sum of the probabilities of all spanning
813 trees of the graph.
815 Parameters
816 ----------
817 G : nx.MultiGraph
818 The undirected support graph for the Held Karp relaxation
820 z : dict
821 The output of `held_karp_ascent()`, a scaled version of the Held-Karp
822 solution.
824 Returns
825 -------
826 gamma : dict
827 The probability distribution which approximately preserves the marginal
828 probabilities of `z`.
829 """
830 from math import exp
831 from math import log as ln
833 def q(e):
834 """
835 The value of q(e) is described in the Asadpour paper is "the
836 probability that edge e will be included in a spanning tree T that is
837 chosen with probability proportional to exp(gamma(T))" which
838 basically means that it is the total probability of the edge appearing
839 across the whole distribution.
841 Parameters
842 ----------
843 e : tuple
844 The `(u, v)` tuple describing the edge we are interested in
846 Returns
847 -------
848 float
849 The probability that a spanning tree chosen according to the
850 current values of gamma will include edge `e`.
851 """
852 # Create the laplacian matrices
853 for u, v, d in G.edges(data=True):
854 d[lambda_key] = exp(gamma[(u, v)])
855 G_Kirchhoff = nx.total_spanning_tree_weight(G, lambda_key)
856 G_e = nx.contracted_edge(G, e, self_loops=False)
857 G_e_Kirchhoff = nx.total_spanning_tree_weight(G_e, lambda_key)
859 # Multiply by the weight of the contracted edge since it is not included
860 # in the total weight of the contracted graph.
861 return exp(gamma[(e[0], e[1])]) * G_e_Kirchhoff / G_Kirchhoff
863 # initialize gamma to the zero dict
864 gamma = {}
865 for u, v, _ in G.edges:
866 gamma[(u, v)] = 0
868 # set epsilon
869 EPSILON = 0.2
871 # pick an edge attribute name that is unlikely to be in the graph
872 lambda_key = "spanning_tree_distribution's secret attribute name for lambda"
874 while True:
875 # We need to know that know that no values of q_e are greater than
876 # (1 + epsilon) * z_e, however changing one gamma value can increase the
877 # value of a different q_e, so we have to complete the for loop without
878 # changing anything for the condition to be meet
879 in_range_count = 0
880 # Search for an edge with q_e > (1 + epsilon) * z_e
881 for u, v in gamma:
882 e = (u, v)
883 q_e = q(e)
884 z_e = z[e]
885 if q_e > (1 + EPSILON) * z_e:
886 delta = ln(
887 (q_e * (1 - (1 + EPSILON / 2) * z_e))
888 / ((1 - q_e) * (1 + EPSILON / 2) * z_e)
889 )
890 gamma[e] -= delta
891 # Check that delta had the desired effect
892 new_q_e = q(e)
893 desired_q_e = (1 + EPSILON / 2) * z_e
894 if round(new_q_e, 8) != round(desired_q_e, 8):
895 raise nx.NetworkXError(
896 f"Unable to modify probability for edge ({u}, {v})"
897 )
898 else:
899 in_range_count += 1
900 # Check if the for loop terminated without changing any gamma
901 if in_range_count == len(gamma):
902 break
904 # Remove the new edge attributes
905 for _, _, d in G.edges(data=True):
906 if lambda_key in d:
907 del d[lambda_key]
909 return gamma
912@nx._dispatch(edge_attrs="weight")
913def greedy_tsp(G, weight="weight", source=None):
914 """Return a low cost cycle starting at `source` and its cost.
916 This approximates a solution to the traveling salesman problem.
917 It finds a cycle of all the nodes that a salesman can visit in order
918 to visit many nodes while minimizing total distance.
919 It uses a simple greedy algorithm.
920 In essence, this function returns a large cycle given a source point
921 for which the total cost of the cycle is minimized.
923 Parameters
924 ----------
925 G : Graph
926 The Graph should be a complete weighted undirected graph.
927 The distance between all pairs of nodes should be included.
929 weight : string, optional (default="weight")
930 Edge data key corresponding to the edge weight.
931 If any edge does not have this attribute the weight is set to 1.
933 source : node, optional (default: first node in list(G))
934 Starting node. If None, defaults to ``next(iter(G))``
936 Returns
937 -------
938 cycle : list of nodes
939 Returns the cycle (list of nodes) that a salesman
940 can follow to minimize total weight of the trip.
942 Raises
943 ------
944 NetworkXError
945 If `G` is not complete, the algorithm raises an exception.
947 Examples
948 --------
949 >>> from networkx.algorithms import approximation as approx
950 >>> G = nx.DiGraph()
951 >>> G.add_weighted_edges_from({
952 ... ("A", "B", 3), ("A", "C", 17), ("A", "D", 14), ("B", "A", 3),
953 ... ("B", "C", 12), ("B", "D", 16), ("C", "A", 13),("C", "B", 12),
954 ... ("C", "D", 4), ("D", "A", 14), ("D", "B", 15), ("D", "C", 2)
955 ... })
956 >>> cycle = approx.greedy_tsp(G, source="D")
957 >>> cost = sum(G[n][nbr]["weight"] for n, nbr in nx.utils.pairwise(cycle))
958 >>> cycle
959 ['D', 'C', 'B', 'A', 'D']
960 >>> cost
961 31
963 Notes
964 -----
965 This implementation of a greedy algorithm is based on the following:
967 - The algorithm adds a node to the solution at every iteration.
968 - The algorithm selects a node not already in the cycle whose connection
969 to the previous node adds the least cost to the cycle.
971 A greedy algorithm does not always give the best solution.
972 However, it can construct a first feasible solution which can
973 be passed as a parameter to an iterative improvement algorithm such
974 as Simulated Annealing, or Threshold Accepting.
976 Time complexity: It has a running time $O(|V|^2)$
977 """
978 # Check that G is a complete graph
979 N = len(G) - 1
980 # This check ignores selfloops which is what we want here.
981 if any(len(nbrdict) - (n in nbrdict) != N for n, nbrdict in G.adj.items()):
982 raise nx.NetworkXError("G must be a complete graph.")
984 if source is None:
985 source = nx.utils.arbitrary_element(G)
987 if G.number_of_nodes() == 2:
988 neighbor = next(G.neighbors(source))
989 return [source, neighbor, source]
991 nodeset = set(G)
992 nodeset.remove(source)
993 cycle = [source]
994 next_node = source
995 while nodeset:
996 nbrdict = G[next_node]
997 next_node = min(nodeset, key=lambda n: nbrdict[n].get(weight, 1))
998 cycle.append(next_node)
999 nodeset.remove(next_node)
1000 cycle.append(cycle[0])
1001 return cycle
1004@py_random_state(9)
1005@nx._dispatch(edge_attrs="weight")
1006def simulated_annealing_tsp(
1007 G,
1008 init_cycle,
1009 weight="weight",
1010 source=None,
1011 temp=100,
1012 move="1-1",
1013 max_iterations=10,
1014 N_inner=100,
1015 alpha=0.01,
1016 seed=None,
1017):
1018 """Returns an approximate solution to the traveling salesman problem.
1020 This function uses simulated annealing to approximate the minimal cost
1021 cycle through the nodes. Starting from a suboptimal solution, simulated
1022 annealing perturbs that solution, occasionally accepting changes that make
1023 the solution worse to escape from a locally optimal solution. The chance
1024 of accepting such changes decreases over the iterations to encourage
1025 an optimal result. In summary, the function returns a cycle starting
1026 at `source` for which the total cost is minimized. It also returns the cost.
1028 The chance of accepting a proposed change is related to a parameter called
1029 the temperature (annealing has a physical analogue of steel hardening
1030 as it cools). As the temperature is reduced, the chance of moves that
1031 increase cost goes down.
1033 Parameters
1034 ----------
1035 G : Graph
1036 `G` should be a complete weighted graph.
1037 The distance between all pairs of nodes should be included.
1039 init_cycle : list of all nodes or "greedy"
1040 The initial solution (a cycle through all nodes returning to the start).
1041 This argument has no default to make you think about it.
1042 If "greedy", use `greedy_tsp(G, weight)`.
1043 Other common starting cycles are `list(G) + [next(iter(G))]` or the final
1044 result of `simulated_annealing_tsp` when doing `threshold_accepting_tsp`.
1046 weight : string, optional (default="weight")
1047 Edge data key corresponding to the edge weight.
1048 If any edge does not have this attribute the weight is set to 1.
1050 source : node, optional (default: first node in list(G))
1051 Starting node. If None, defaults to ``next(iter(G))``
1053 temp : int, optional (default=100)
1054 The algorithm's temperature parameter. It represents the initial
1055 value of temperature
1057 move : "1-1" or "1-0" or function, optional (default="1-1")
1058 Indicator of what move to use when finding new trial solutions.
1059 Strings indicate two special built-in moves:
1061 - "1-1": 1-1 exchange which transposes the position
1062 of two elements of the current solution.
1063 The function called is :func:`swap_two_nodes`.
1064 For example if we apply 1-1 exchange in the solution
1065 ``A = [3, 2, 1, 4, 3]``
1066 we can get the following by the transposition of 1 and 4 elements:
1067 ``A' = [3, 2, 4, 1, 3]``
1068 - "1-0": 1-0 exchange which moves an node in the solution
1069 to a new position.
1070 The function called is :func:`move_one_node`.
1071 For example if we apply 1-0 exchange in the solution
1072 ``A = [3, 2, 1, 4, 3]``
1073 we can transfer the fourth element to the second position:
1074 ``A' = [3, 4, 2, 1, 3]``
1076 You may provide your own functions to enact a move from
1077 one solution to a neighbor solution. The function must take
1078 the solution as input along with a `seed` input to control
1079 random number generation (see the `seed` input here).
1080 Your function should maintain the solution as a cycle with
1081 equal first and last node and all others appearing once.
1082 Your function should return the new solution.
1084 max_iterations : int, optional (default=10)
1085 Declared done when this number of consecutive iterations of
1086 the outer loop occurs without any change in the best cost solution.
1088 N_inner : int, optional (default=100)
1089 The number of iterations of the inner loop.
1091 alpha : float between (0, 1), optional (default=0.01)
1092 Percentage of temperature decrease in each iteration
1093 of outer loop
1095 seed : integer, random_state, or None (default)
1096 Indicator of random number generation state.
1097 See :ref:`Randomness<randomness>`.
1099 Returns
1100 -------
1101 cycle : list of nodes
1102 Returns the cycle (list of nodes) that a salesman
1103 can follow to minimize total weight of the trip.
1105 Raises
1106 ------
1107 NetworkXError
1108 If `G` is not complete the algorithm raises an exception.
1110 Examples
1111 --------
1112 >>> from networkx.algorithms import approximation as approx
1113 >>> G = nx.DiGraph()
1114 >>> G.add_weighted_edges_from({
1115 ... ("A", "B", 3), ("A", "C", 17), ("A", "D", 14), ("B", "A", 3),
1116 ... ("B", "C", 12), ("B", "D", 16), ("C", "A", 13),("C", "B", 12),
1117 ... ("C", "D", 4), ("D", "A", 14), ("D", "B", 15), ("D", "C", 2)
1118 ... })
1119 >>> cycle = approx.simulated_annealing_tsp(G, "greedy", source="D")
1120 >>> cost = sum(G[n][nbr]["weight"] for n, nbr in nx.utils.pairwise(cycle))
1121 >>> cycle
1122 ['D', 'C', 'B', 'A', 'D']
1123 >>> cost
1124 31
1125 >>> incycle = ["D", "B", "A", "C", "D"]
1126 >>> cycle = approx.simulated_annealing_tsp(G, incycle, source="D")
1127 >>> cost = sum(G[n][nbr]["weight"] for n, nbr in nx.utils.pairwise(cycle))
1128 >>> cycle
1129 ['D', 'C', 'B', 'A', 'D']
1130 >>> cost
1131 31
1133 Notes
1134 -----
1135 Simulated Annealing is a metaheuristic local search algorithm.
1136 The main characteristic of this algorithm is that it accepts
1137 even solutions which lead to the increase of the cost in order
1138 to escape from low quality local optimal solutions.
1140 This algorithm needs an initial solution. If not provided, it is
1141 constructed by a simple greedy algorithm. At every iteration, the
1142 algorithm selects thoughtfully a neighbor solution.
1143 Consider $c(x)$ cost of current solution and $c(x')$ cost of a
1144 neighbor solution.
1145 If $c(x') - c(x) <= 0$ then the neighbor solution becomes the current
1146 solution for the next iteration. Otherwise, the algorithm accepts
1147 the neighbor solution with probability $p = exp - ([c(x') - c(x)] / temp)$.
1148 Otherwise the current solution is retained.
1150 `temp` is a parameter of the algorithm and represents temperature.
1152 Time complexity:
1153 For $N_i$ iterations of the inner loop and $N_o$ iterations of the
1154 outer loop, this algorithm has running time $O(N_i * N_o * |V|)$.
1156 For more information and how the algorithm is inspired see:
1157 http://en.wikipedia.org/wiki/Simulated_annealing
1158 """
1159 if move == "1-1":
1160 move = swap_two_nodes
1161 elif move == "1-0":
1162 move = move_one_node
1163 if init_cycle == "greedy":
1164 # Construct an initial solution using a greedy algorithm.
1165 cycle = greedy_tsp(G, weight=weight, source=source)
1166 if G.number_of_nodes() == 2:
1167 return cycle
1169 else:
1170 cycle = list(init_cycle)
1171 if source is None:
1172 source = cycle[0]
1173 elif source != cycle[0]:
1174 raise nx.NetworkXError("source must be first node in init_cycle")
1175 if cycle[0] != cycle[-1]:
1176 raise nx.NetworkXError("init_cycle must be a cycle. (return to start)")
1178 if len(cycle) - 1 != len(G) or len(set(G.nbunch_iter(cycle))) != len(G):
1179 raise nx.NetworkXError("init_cycle should be a cycle over all nodes in G.")
1181 # Check that G is a complete graph
1182 N = len(G) - 1
1183 # This check ignores selfloops which is what we want here.
1184 if any(len(nbrdict) - (n in nbrdict) != N for n, nbrdict in G.adj.items()):
1185 raise nx.NetworkXError("G must be a complete graph.")
1187 if G.number_of_nodes() == 2:
1188 neighbor = next(G.neighbors(source))
1189 return [source, neighbor, source]
1191 # Find the cost of initial solution
1192 cost = sum(G[u][v].get(weight, 1) for u, v in pairwise(cycle))
1194 count = 0
1195 best_cycle = cycle.copy()
1196 best_cost = cost
1197 while count <= max_iterations and temp > 0:
1198 count += 1
1199 for i in range(N_inner):
1200 adj_sol = move(cycle, seed)
1201 adj_cost = sum(G[u][v].get(weight, 1) for u, v in pairwise(adj_sol))
1202 delta = adj_cost - cost
1203 if delta <= 0:
1204 # Set current solution the adjacent solution.
1205 cycle = adj_sol
1206 cost = adj_cost
1208 if cost < best_cost:
1209 count = 0
1210 best_cycle = cycle.copy()
1211 best_cost = cost
1212 else:
1213 # Accept even a worse solution with probability p.
1214 p = math.exp(-delta / temp)
1215 if p >= seed.random():
1216 cycle = adj_sol
1217 cost = adj_cost
1218 temp -= temp * alpha
1220 return best_cycle
1223@py_random_state(9)
1224@nx._dispatch(edge_attrs="weight")
1225def threshold_accepting_tsp(
1226 G,
1227 init_cycle,
1228 weight="weight",
1229 source=None,
1230 threshold=1,
1231 move="1-1",
1232 max_iterations=10,
1233 N_inner=100,
1234 alpha=0.1,
1235 seed=None,
1236):
1237 """Returns an approximate solution to the traveling salesman problem.
1239 This function uses threshold accepting methods to approximate the minimal cost
1240 cycle through the nodes. Starting from a suboptimal solution, threshold
1241 accepting methods perturb that solution, accepting any changes that make
1242 the solution no worse than increasing by a threshold amount. Improvements
1243 in cost are accepted, but so are changes leading to small increases in cost.
1244 This allows the solution to leave suboptimal local minima in solution space.
1245 The threshold is decreased slowly as iterations proceed helping to ensure
1246 an optimum. In summary, the function returns a cycle starting at `source`
1247 for which the total cost is minimized.
1249 Parameters
1250 ----------
1251 G : Graph
1252 `G` should be a complete weighted graph.
1253 The distance between all pairs of nodes should be included.
1255 init_cycle : list or "greedy"
1256 The initial solution (a cycle through all nodes returning to the start).
1257 This argument has no default to make you think about it.
1258 If "greedy", use `greedy_tsp(G, weight)`.
1259 Other common starting cycles are `list(G) + [next(iter(G))]` or the final
1260 result of `simulated_annealing_tsp` when doing `threshold_accepting_tsp`.
1262 weight : string, optional (default="weight")
1263 Edge data key corresponding to the edge weight.
1264 If any edge does not have this attribute the weight is set to 1.
1266 source : node, optional (default: first node in list(G))
1267 Starting node. If None, defaults to ``next(iter(G))``
1269 threshold : int, optional (default=1)
1270 The algorithm's threshold parameter. It represents the initial
1271 threshold's value
1273 move : "1-1" or "1-0" or function, optional (default="1-1")
1274 Indicator of what move to use when finding new trial solutions.
1275 Strings indicate two special built-in moves:
1277 - "1-1": 1-1 exchange which transposes the position
1278 of two elements of the current solution.
1279 The function called is :func:`swap_two_nodes`.
1280 For example if we apply 1-1 exchange in the solution
1281 ``A = [3, 2, 1, 4, 3]``
1282 we can get the following by the transposition of 1 and 4 elements:
1283 ``A' = [3, 2, 4, 1, 3]``
1284 - "1-0": 1-0 exchange which moves an node in the solution
1285 to a new position.
1286 The function called is :func:`move_one_node`.
1287 For example if we apply 1-0 exchange in the solution
1288 ``A = [3, 2, 1, 4, 3]``
1289 we can transfer the fourth element to the second position:
1290 ``A' = [3, 4, 2, 1, 3]``
1292 You may provide your own functions to enact a move from
1293 one solution to a neighbor solution. The function must take
1294 the solution as input along with a `seed` input to control
1295 random number generation (see the `seed` input here).
1296 Your function should maintain the solution as a cycle with
1297 equal first and last node and all others appearing once.
1298 Your function should return the new solution.
1300 max_iterations : int, optional (default=10)
1301 Declared done when this number of consecutive iterations of
1302 the outer loop occurs without any change in the best cost solution.
1304 N_inner : int, optional (default=100)
1305 The number of iterations of the inner loop.
1307 alpha : float between (0, 1), optional (default=0.1)
1308 Percentage of threshold decrease when there is at
1309 least one acceptance of a neighbor solution.
1310 If no inner loop moves are accepted the threshold remains unchanged.
1312 seed : integer, random_state, or None (default)
1313 Indicator of random number generation state.
1314 See :ref:`Randomness<randomness>`.
1316 Returns
1317 -------
1318 cycle : list of nodes
1319 Returns the cycle (list of nodes) that a salesman
1320 can follow to minimize total weight of the trip.
1322 Raises
1323 ------
1324 NetworkXError
1325 If `G` is not complete the algorithm raises an exception.
1327 Examples
1328 --------
1329 >>> from networkx.algorithms import approximation as approx
1330 >>> G = nx.DiGraph()
1331 >>> G.add_weighted_edges_from({
1332 ... ("A", "B", 3), ("A", "C", 17), ("A", "D", 14), ("B", "A", 3),
1333 ... ("B", "C", 12), ("B", "D", 16), ("C", "A", 13),("C", "B", 12),
1334 ... ("C", "D", 4), ("D", "A", 14), ("D", "B", 15), ("D", "C", 2)
1335 ... })
1336 >>> cycle = approx.threshold_accepting_tsp(G, "greedy", source="D")
1337 >>> cost = sum(G[n][nbr]["weight"] for n, nbr in nx.utils.pairwise(cycle))
1338 >>> cycle
1339 ['D', 'C', 'B', 'A', 'D']
1340 >>> cost
1341 31
1342 >>> incycle = ["D", "B", "A", "C", "D"]
1343 >>> cycle = approx.threshold_accepting_tsp(G, incycle, source="D")
1344 >>> cost = sum(G[n][nbr]["weight"] for n, nbr in nx.utils.pairwise(cycle))
1345 >>> cycle
1346 ['D', 'C', 'B', 'A', 'D']
1347 >>> cost
1348 31
1350 Notes
1351 -----
1352 Threshold Accepting is a metaheuristic local search algorithm.
1353 The main characteristic of this algorithm is that it accepts
1354 even solutions which lead to the increase of the cost in order
1355 to escape from low quality local optimal solutions.
1357 This algorithm needs an initial solution. This solution can be
1358 constructed by a simple greedy algorithm. At every iteration, it
1359 selects thoughtfully a neighbor solution.
1360 Consider $c(x)$ cost of current solution and $c(x')$ cost of
1361 neighbor solution.
1362 If $c(x') - c(x) <= threshold$ then the neighbor solution becomes the current
1363 solution for the next iteration, where the threshold is named threshold.
1365 In comparison to the Simulated Annealing algorithm, the Threshold
1366 Accepting algorithm does not accept very low quality solutions
1367 (due to the presence of the threshold value). In the case of
1368 Simulated Annealing, even a very low quality solution can
1369 be accepted with probability $p$.
1371 Time complexity:
1372 It has a running time $O(m * n * |V|)$ where $m$ and $n$ are the number
1373 of times the outer and inner loop run respectively.
1375 For more information and how algorithm is inspired see:
1376 https://doi.org/10.1016/0021-9991(90)90201-B
1378 See Also
1379 --------
1380 simulated_annealing_tsp
1382 """
1383 if move == "1-1":
1384 move = swap_two_nodes
1385 elif move == "1-0":
1386 move = move_one_node
1387 if init_cycle == "greedy":
1388 # Construct an initial solution using a greedy algorithm.
1389 cycle = greedy_tsp(G, weight=weight, source=source)
1390 if G.number_of_nodes() == 2:
1391 return cycle
1393 else:
1394 cycle = list(init_cycle)
1395 if source is None:
1396 source = cycle[0]
1397 elif source != cycle[0]:
1398 raise nx.NetworkXError("source must be first node in init_cycle")
1399 if cycle[0] != cycle[-1]:
1400 raise nx.NetworkXError("init_cycle must be a cycle. (return to start)")
1402 if len(cycle) - 1 != len(G) or len(set(G.nbunch_iter(cycle))) != len(G):
1403 raise nx.NetworkXError("init_cycle is not all and only nodes.")
1405 # Check that G is a complete graph
1406 N = len(G) - 1
1407 # This check ignores selfloops which is what we want here.
1408 if any(len(nbrdict) - (n in nbrdict) != N for n, nbrdict in G.adj.items()):
1409 raise nx.NetworkXError("G must be a complete graph.")
1411 if G.number_of_nodes() == 2:
1412 neighbor = list(G.neighbors(source))[0]
1413 return [source, neighbor, source]
1415 # Find the cost of initial solution
1416 cost = sum(G[u][v].get(weight, 1) for u, v in pairwise(cycle))
1418 count = 0
1419 best_cycle = cycle.copy()
1420 best_cost = cost
1421 while count <= max_iterations:
1422 count += 1
1423 accepted = False
1424 for i in range(N_inner):
1425 adj_sol = move(cycle, seed)
1426 adj_cost = sum(G[u][v].get(weight, 1) for u, v in pairwise(adj_sol))
1427 delta = adj_cost - cost
1428 if delta <= threshold:
1429 accepted = True
1431 # Set current solution the adjacent solution.
1432 cycle = adj_sol
1433 cost = adj_cost
1435 if cost < best_cost:
1436 count = 0
1437 best_cycle = cycle.copy()
1438 best_cost = cost
1439 if accepted:
1440 threshold -= threshold * alpha
1442 return best_cycle