"""
=================================
旅行商问题 (TSP)
=================================
实现近似算法来解决和近似TSP问题。
已实现的算法类别包括：
- 克里斯托费德斯算法（提供3/2近似解）
- 贪心算法
- 模拟退火算法 (SA)
- 阈值接受算法 (TA)
- Asadpour 非对称旅行商算法
旅行商问题试图找到，在给定销售员必须访问的所有点之间的权重（距离）的情况下，满足以下条件的路线：
- 销售员旅行的总距离（成本）最小化。
- 销售员返回起点。
- 注意，对于完全图，销售员每个点只访问一次。
函数 `travelling_salesman_problem` 通过找到所有点对之间的最短路径，有效地将问题转换为完全图问题，从而允许处理不完全图。它在该问题上调用一种近似方法，然后使用之前找到的最短路径将结果转换回原始图。
TSP 是组合优化中的一个 NP 难问题，在运筹学和理论计算机科学中非常重要。
http://en.wikipedia.org/wiki/Travelling_salesman_problem
"""
import math
import networkx as nx
from networkx.algorithms.tree.mst import random_spanning_tree
from networkx.utils import not_implemented_for, pairwise, py_random_state
__all__ = [
"traveling_salesman_problem",
"christofides",
"asadpour_atsp",
"greedy_tsp",
"simulated_annealing_tsp",
"threshold_accepting_tsp",
]
def swap_two_nodes(soln, seed):
"""交换 `soln` 中的两个节点以得到一个邻域解。
Parameters
----------
soln : 节点列表
当前的节点循环
seed : 整数, random_state, 或 None (默认)
随机数生成状态的指示器。
参见 :ref:`随机性<randomness>` 。
Returns
-------
列表
应用移动后的解。（一个邻域解。）
Notes
-----
此函数假设传入的列表 `soln` 是一个循环
（即第一个和最后一个元素相同），并且
我们不希望任何移动改变列表中的第一个节点
（因此也不会改变最后一个节点）。
输入列表也会被改变并返回。如果需要，请复制一份。
See Also
--------
move_one_node
"""
a, b = seed.sample(range(1, len(soln) - 1), k=2)
soln[a], soln[b] = soln[b], soln[a]
return soln
def move_one_node(soln, seed):
"""将一个节点移动到另一个位置以生成一个邻域解。
要移动的节点和目标位置是随机选择的。
第一个和最后一个节点保持不变，因为解必须是从该节点开始的循环。
Parameters
----------
soln : 节点列表
当前的节点循环
seed : 整数, random_state, 或 None (默认)
随机数生成状态的指示器。
参见 :ref:`随机性<randomness>` 。
Returns
-------
列表
应用移动后的解。（一个邻域解。）
Notes
-----
此函数假设传入的列表 `soln` 是一个循环
（即第一个和最后一个元素相同），并且我们不希望任何移动改变列表中的第一个节点
（因此也不改变最后一个节点）。
输入列表会被修改并返回。如果需要，请复制一份。
See Also
--------
swap_two_nodes
"""
a, b = seed.sample(range(1, len(soln) - 1), k=2)
soln.insert(b, soln.pop(a))
return soln
[docs]
@not_implemented_for("directed")
@nx._dispatchable(edge_attrs="weight")
def christofides(G, weight="weight", tree=None):
"""近似求解旅行商问题
计算完全无向图中旅行商问题的3/2近似解
使用Christofides [1]_算法。
Parameters
----------
G : 图
`G` 应为一个完全加权无向图。
应包含所有节点对之间的距离。
weight : 字符串, 可选 (默认="weight")
对应于边权的边数据键。
如果任何边没有此属性，则权重设置为1。
tree : NetworkX图或None (默认: None)
图G的最小生成树。或者，如果为None，则使用
:func:`networkx.minimum_spanning_tree` 计算最小生成树。
Returns
-------
列表
`G` 中沿着一个3/2近似最小哈密顿环的节点列表。
References
----------
.. [1] Christofides, Nicos. "旅行商问题的新启发式算法的
最坏情况分析." No. RR-388. Carnegie-Mellon Univ
Pittsburgh Pa 管理科学研究组, 1976.
"""
# Remove selfloops if necessary
loop_nodes = nx.nodes_with_selfloops(G)
try:
node = next(loop_nodes)
except StopIteration:
pass
else:
G = G.copy()
G.remove_edge(node, node)
G.remove_edges_from((n, n) for n in loop_nodes)
# Check that G is a complete graph
N = len(G) - 1
# This check ignores selfloops which is what we want here.
if any(len(nbrdict) != N for n, nbrdict in G.adj.items()):
raise nx.NetworkXError("G must be a complete graph.")
if tree is None:
tree = nx.minimum_spanning_tree(G, weight=weight)
L = G.copy()
L.remove_nodes_from([v for v, degree in tree.degree if not (degree % 2)])
MG = nx.MultiGraph()
MG.add_edges_from(tree.edges)
edges = nx.min_weight_matching(L, weight=weight)
MG.add_edges_from(edges)
return _shortcutting(nx.eulerian_circuit(MG))
def _shortcutting(circuit):
"""移除路径中的重复节点"""
nodes = []
for u, v in circuit:
if v in nodes:
continue
if not nodes:
nodes.append(u)
nodes.append(v)
nodes.append(nodes[0])
return nodes
[docs]
@nx._dispatchable(edge_attrs="weight")
def traveling_salesman_problem(
G, weight="weight", nodes=None, cycle=True, method=None, **kwargs
):
"""在 `G` 中找到连接指定节点的最短路径
此函数允许在非完全图和/或不需要访问所有节点的网络中近似解决旅行商问题。
此函数分两步进行。首先，它使用 `nodes` 中节点之间的所有对最短路径创建一个完全图。新图中的边权重是原始图中每对节点之间路径的长度。
其次，使用一种算法（默认：无向图使用 `christofides` ，有向图使用 `asadpour_atsp` ）来近似这个新图上的最小哈密顿循环。可用的算法有：
- christofides
- greedy_tsp
- simulated_annealing_tsp
- threshold_accepting_tsp
- asadpour_atsp
一旦找到哈密顿循环，此函数会进行后处理以适应原始图的结构。如果 `cycle` 为 ``False`` ，则移除最大权重边以形成哈密顿路径。然后，用于分析的新完全图中的每条边都被替换为原始图中这些节点之间的最短路径。如果输入图 `G` 包含不遵循三角不等式的边权重，例如当 `G` 不是完全图时（即不存在的边的长度为无穷大），则返回的路径可能包含一些重复节点（起始节点除外）。
Parameters
----------
G : NetworkX 图
一个可能有权重的图
nodes : 节点集合（默认=G.nodes）
要访问的节点集合（列表、集合等）
weight : 字符串，可选（默认="weight"）
对应边权重的边数据键。如果任何边没有此属性，则权重设置为1。
cycle : 布尔值（默认: True）
指示是否应返回一个循环，还是一个路径。注意：该循环是近似的最小循环。路径只是移除该循环中的最大边。
method : 函数（默认: None）
一个返回所有节点上循环并近似解决完全图上旅行商问题的函数。返回的循环用于在 `G` 上找到相应的解决方案。 `method` 应该是可调用的；接受输入 `G` 和 `weight` ；并返回沿循环的节点列表。
提供的选项包括 :func:`christofides` , :func:`greedy_tsp` , :func:`simulated_annealing_tsp` 和 :func:`threshold_accepting_tsp` 。
如果 `method` 为 None：对于无向图 `G` 使用 :func:`christofides` ，对于有向图 `G` 使用 :func:`asadpour_atsp` 。
**kwargs : 字典
传递给传入的 `method` 函数的其他关键字参数。
Returns
-------
列表
`G` 中沿近似最小路径通过 `nodes` 的节点列表。
Raises
------
NetworkXError
如果 `G` 是有向图，它必须是强连通的，否则无法生成完全图版本。
Examples
--------
>>> tsp = nx.approximation.traveling_salesman_problem
>>> G = nx.cycle_graph(9)
>>> G[4][5]["weight"] = 5 # 所有其他权重为1
>>> tsp(G, nodes=[3, 6])
[3, 2, 1, 0, 8, 7, 6, 7, 8, 0, 1, 2, 3]
>>> path = tsp(G, cycle=False)
>>> path in ([4, 3, 2, 1, 0, 8, 7, 6, 5], [5, 6, 7, 8, 0, 1, 2, 3, 4])
True
虽然不再需要，你仍然可以构建（柯里化）你自己的函数来为方法提供参数值。
>>> SA_tsp = nx.approximation.simulated_annealing_tsp
>>> method = lambda G, weight: SA_tsp(G, "greedy", weight=weight, temp=500)
>>> path = tsp(G, cycle=False, method=method)
>>> path in ([4, 3, 2, 1, 0, 8, 7, 6, 5], [5, 6, 7, 8, 0, 1, 2, 3, 4])
True
否则，直接将其他关键字参数传递给 tsp 函数。
>>> path = tsp(
... G,
... cycle=False,
... method=nx.approximation.simulated_annealing_tsp,
... init_cycle="greedy",
... temp=500,
... )
>>> path in ([4, 3, 2, 1, 0, 8, 7, 6, 5], [5, 6, 7, 8, 0, 1, 2, 3, 4])
True
"""
if method is None:
if G.is_directed():
method = asadpour_atsp
else:
method = christofides
if nodes is None:
nodes = list(G.nodes)
dist = {}
path = {}
for n, (d, p) in nx.all_pairs_dijkstra(G, weight=weight):
dist[n] = d
path[n] = p
if G.is_directed():
# If the graph is not strongly connected, raise an exception
if not nx.is_strongly_connected(G):
raise nx.NetworkXError("G is not strongly connected")
GG = nx.DiGraph()
else:
GG = nx.Graph()
for u in nodes:
for v in nodes:
if u == v:
continue
GG.add_edge(u, v, weight=dist[u][v])
best_GG = method(GG, weight=weight, **kwargs)
if not cycle:
# find and remove the biggest edge
(u, v) = max(pairwise(best_GG), key=lambda x: dist[x[0]][x[1]])
pos = best_GG.index(u) + 1
while best_GG[pos] != v:
pos = best_GG[pos:].index(u) + 1
best_GG = best_GG[pos:-1] + best_GG[:pos]
best_path = []
for u, v in pairwise(best_GG):
best_path.extend(path[u][v][:-1])
best_path.append(v)
return best_path
[docs]
@not_implemented_for("undirected")
@py_random_state(2)
@nx._dispatchable(edge_attrs="weight", mutates_input=True)
def asadpour_atsp(G, weight="weight", seed=None, source=None):
"""返回旅行商问题的近似解。
这个近似解是Asadpour等人开发的非对称旅行商问题中已知的最佳近似算法之一，[1]。该算法首先解决Held-Karp松弛问题，以找到循环权重的下界。接下来，它构建一个无向生成树的指数分布，其中树中边的概率与该边的权重相对应，使用最大熵舍入方案。然后我们采样该分布$2 \\lceil \\ln n \\rceil$次，并在弧的方向添加回边后保存最小采样树。最后，我们增强并短路该图，以找到销售员的近似旅行路线。
Parameters
----------
G : nx.DiGraph
图应为一个完整的加权有向图。所有节点对之间的距离应包括在内，并且应满足三角不等式。也就是说，任意两个节点之间的直接边应为成本最低的路径。
weight : string, 可选 (默认="weight")
对应于边权重的边数据键。
如果任何边没有此属性，则权重设置为1。
seed : integer, random_state, 或 None (默认)
随机数生成状态的指示符。
参见 :ref:`Randomness<randomness>`。
source : 节点标签 (默认=`None`)
如果给定，返回从给定节点开始和结束的循环。
Returns
-------
cycle : 节点列表
返回销售员可以遵循以最小化旅行总权重的循环（节点列表）。
Raises
------
NetworkXError
如果`G`不完整或节点数少于两个，算法会引发异常。
NetworkXError
如果`source`不是`None`且不是`G`中的节点，算法会引发异常。
NetworkXNotImplemented
如果`G`是无向图。
References
----------
.. [1] A. Asadpour, M. X. Goemans, A. Madry, S. O. Gharan, 和 A. Saberi,
非对称旅行商问题的o(log n/log log n)近似算法, 运筹学, 65 (2017),
pp. 1043–1061
Examples
--------
>>> import networkx as nx
>>> import networkx.algorithms.approximation as approx
>>> G = nx.complete_graph(3, create_using=nx.DiGraph)
>>> nx.set_edge_attributes(
... G,
... {(0, 1): 2, (1, 2): 2, (2, 0): 2, (0, 2): 1, (2, 1): 1, (1, 0): 1},
... "weight",
... )
>>> tour = approx.asadpour_atsp(G, source=0)
>>> tour
[0, 2, 1, 0]
"""
from math import ceil, exp
from math import log as ln
# Check that G is a complete graph
N = len(G) - 1
if N < 2:
raise nx.NetworkXError("G must have at least two nodes")
# This check ignores selfloops which is what we want here.
if any(len(nbrdict) - (n in nbrdict) != N for n, nbrdict in G.adj.items()):
raise nx.NetworkXError("G is not a complete DiGraph")
# Check that the source vertex, if given, is in the graph
if source is not None and source not in G.nodes:
raise nx.NetworkXError("Given source node not in G.")
opt_hk, z_star = held_karp_ascent(G, weight)
# Test to see if the ascent method found an integer solution or a fractional
# solution. If it is integral then z_star is a nx.Graph, otherwise it is
# a dict
if not isinstance(z_star, dict):
# Here we are using the shortcutting method to go from the list of edges
# returned from eulerian_circuit to a list of nodes
return _shortcutting(nx.eulerian_circuit(z_star, source=source))
# Create the undirected support of z_star
z_support = nx.MultiGraph()
for u, v in z_star:
if (u, v) not in z_support.edges:
edge_weight = min(G[u][v][weight], G[v][u][weight])
z_support.add_edge(u, v, **{weight: edge_weight})
# Create the exponential distribution of spanning trees
gamma = spanning_tree_distribution(z_support, z_star)
# Write the lambda values to the edges of z_support
z_support = nx.Graph(z_support)
lambda_dict = {(u, v): exp(gamma[(u, v)]) for u, v in z_support.edges()}
nx.set_edge_attributes(z_support, lambda_dict, "weight")
del gamma, lambda_dict
# Sample 2 * ceil( ln(n) ) spanning trees and record the minimum one
minimum_sampled_tree = None
minimum_sampled_tree_weight = math.inf
for _ in range(2 * ceil(ln(G.number_of_nodes()))):
sampled_tree = random_spanning_tree(z_support, "weight", seed=seed)
sampled_tree_weight = sampled_tree.size(weight)
if sampled_tree_weight < minimum_sampled_tree_weight:
minimum_sampled_tree = sampled_tree.copy()
minimum_sampled_tree_weight = sampled_tree_weight
# Orient the edges in that tree to keep the cost of the tree the same.
t_star = nx.MultiDiGraph()
for u, v, d in minimum_sampled_tree.edges(data=weight):
if d == G[u][v][weight]:
t_star.add_edge(u, v, **{weight: d})
else:
t_star.add_edge(v, u, **{weight: d})
# Find the node demands needed to neutralize the flow of t_star in G
node_demands = {n: t_star.out_degree(n) - t_star.in_degree(n) for n in t_star}
nx.set_node_attributes(G, node_demands, "demand")
# Find the min_cost_flow
flow_dict = nx.min_cost_flow(G, "demand")
# Build the flow into t_star
for source, values in flow_dict.items():
for target in values:
if (source, target) not in t_star.edges and values[target] > 0:
# IF values[target] > 0 we have to add that many edges
for _ in range(values[target]):
t_star.add_edge(source, target)
# Return the shortcut eulerian circuit
circuit = nx.eulerian_circuit(t_star, source=source)
return _shortcutting(circuit)
@nx._dispatchable(edge_attrs="weight", mutates_input=True, returns_graph=True)
def held_karp_ascent(G, weight="weight"):
"""最小化TSP的Held-Karp松弛问题
解决输入完全有向图的Held-Karp松弛问题，并缩放输出解以用于Asadpour [1]_ ASTP算法。
Held-Karp松弛定义了ATSP解决方案的下界，尽管它确实返回了一个分数解。这在Asadpour算法中用作初始解，随后在生成树多面体内舍入为整数树。此函数使用[2]_中的分支定界方法解决松弛问题。
Parameters
----------
G : nx.DiGraph
图应为完全加权有向图。
所有节点对之间的距离应包含在内。
weight : string, 可选 (默认="weight")
对应于边权重的边数据键。
如果任何边没有此属性，则权重设置为1。
Returns
-------
OPT : float
Held-Karp松弛最优解的成本
z : dict 或 nx.Graph
Asadpour算法中使用的Held-Karp松弛最优解的对称化和缩放版本。
如果找到整数解，则该解为ATSP问题的最优解，并返回该解。
References
----------
.. [1] A. Asadpour, M. X. Goemans, A. Madry, S. O. Gharan, 和 A. Saberi,
不对称旅行商问题的o(log n/log log n)近似算法, 运筹学, 65 (2017),
pp. 1043–1061
.. [2] M. Held, R. M. Karp, 旅行商问题与最小生成树, 运筹学, 1970-11-01, Vol. 18 (6),
pp.1138-1162
"""
import numpy as np
from scipy import optimize
def k_pi():
"""找到图G在点pi处的最小1-树形结构集合。
Returns
-------
集合
最小1-树形结构集合
"""
# Create a copy of G without vertex 1.
G_1 = G.copy()
minimum_1_arborescences = set()
minimum_1_arborescence_weight = math.inf
# node is node '1' in the Held and Karp paper
n = next(G.__iter__())
G_1.remove_node(n)
# Iterate over the spanning arborescences of the graph until we know
# that we have found the minimum 1-arborescences. My proposed strategy
# is to find the most extensive root to connect to from 'node 1' and
# the least expensive one. We then iterate over arborescences until
# the cost of the basic arborescence is the cost of the minimum one
# plus the difference between the most and least expensive roots,
# that way the cost of connecting 'node 1' will by definition not by
# minimum
min_root = {"node": None, weight: math.inf}
max_root = {"node": None, weight: -math.inf}
for u, v, d in G.edges(n, data=True):
if d[weight] < min_root[weight]:
min_root = {"node": v, weight: d[weight]}
if d[weight] > max_root[weight]:
max_root = {"node": v, weight: d[weight]}
min_in_edge = min(G.in_edges(n, data=True), key=lambda x: x[2][weight])
min_root[weight] = min_root[weight] + min_in_edge[2][weight]
max_root[weight] = max_root[weight] + min_in_edge[2][weight]
min_arb_weight = math.inf
for arb in nx.ArborescenceIterator(G_1):
arb_weight = arb.size(weight)
if min_arb_weight == math.inf:
min_arb_weight = arb_weight
elif arb_weight > min_arb_weight + max_root[weight] - min_root[weight]:
break
# We have to pick the root node of the arborescence for the out
# edge of the first vertex as that is the only node without an
# edge directed into it.
for N, deg in arb.in_degree:
if deg == 0:
# root found
arb.add_edge(n, N, **{weight: G[n][N][weight]})
arb_weight += G[n][N][weight]
break
# We can pick the minimum weight in-edge for the vertex with
# a cycle. If there are multiple edges with the same, minimum
# weight, We need to add all of them.
#
# Delete the edge (N, v) so that we cannot pick it.
edge_data = G[N][n]
G.remove_edge(N, n)
min_weight = min(G.in_edges(n, data=weight), key=lambda x: x[2])[2]
min_edges = [
(u, v, d) for u, v, d in G.in_edges(n, data=weight) if d == min_weight
]
for u, v, d in min_edges:
new_arb = arb.copy()
new_arb.add_edge(u, v, **{weight: d})
new_arb_weight = arb_weight + d
# Check to see the weight of the arborescence, if it is a
# new minimum, clear all of the old potential minimum
# 1-arborescences and add this is the only one. If its
# weight is above the known minimum, do not add it.
if new_arb_weight < minimum_1_arborescence_weight:
minimum_1_arborescences.clear()
minimum_1_arborescence_weight = new_arb_weight
# We have a 1-arborescence, add it to the set
if new_arb_weight == minimum_1_arborescence_weight:
minimum_1_arborescences.add(new_arb)
G.add_edge(N, n, **edge_data)
return minimum_1_arborescences
def direction_of_ascent():
"""在点 pi 处找到上升方向。
更多信息请参见 [1]_。
Returns
-------
dict
从图的节点映射，表示上升方向。
References
----------
.. [1] M. Held, R. M. Karp, 旅行推销员问题与最小生成树，运筹学，1970-11-01，第18卷（6），
第1138-1162页
"""
# 1. Set d equal to the zero n-vector.
d = {}
for n in G:
d[n] = 0
del n
# 2. Find a 1-Arborescence T^k such that k is in K(pi, d).
minimum_1_arborescences = k_pi()
while True:
# Reduce K(pi) to K(pi, d)
# Find the arborescence in K(pi) which increases the lest in
# direction d
min_k_d_weight = math.inf
min_k_d = None
for arborescence in minimum_1_arborescences:
weighted_cost = 0
for n, deg in arborescence.degree:
weighted_cost += d[n] * (deg - 2)
if weighted_cost < min_k_d_weight:
min_k_d_weight = weighted_cost
min_k_d = arborescence
# 3. If sum of d_i * v_{i, k} is greater than zero, terminate
if min_k_d_weight > 0:
return d, min_k_d
# 4. d_i = d_i + v_{i, k}
for n, deg in min_k_d.degree:
d[n] += deg - 2
# Check that we do not need to terminate because the direction
# of ascent does not exist. This is done with linear
# programming.
c = np.full(len(minimum_1_arborescences), -1, dtype=int)
a_eq = np.empty((len(G) + 1, len(minimum_1_arborescences)), dtype=int)
b_eq = np.zeros(len(G) + 1, dtype=int)
b_eq[len(G)] = 1
for arb_count, arborescence in enumerate(minimum_1_arborescences):
n_count = len(G) - 1
for n, deg in arborescence.degree:
a_eq[n_count][arb_count] = deg - 2
n_count -= 1
a_eq[len(G)][arb_count] = 1
program_result = optimize.linprog(
c, A_eq=a_eq, b_eq=b_eq, method="highs-ipm"
)
# If the constants exist, then the direction of ascent doesn't
if program_result.success:
# There is no direction of ascent
return None, minimum_1_arborescences
# 5. GO TO 2
def find_epsilon(k, d):
"""给定在 pi 处的上升方向，找到在该方向上我们可以行进的最大距离。
Parameters
----------
k_xy : set
具有在上升方向上最小增加率的 1-树形集合
d : dict
上升方向
Returns
-------
float
在方向 `d` 上我们可以行进的距离
"""
min_epsilon = math.inf
for e_u, e_v, e_w in G.edges(data=weight):
if (e_u, e_v) in k.edges:
continue
# Now, I have found a condition which MUST be true for the edges to
# be a valid substitute. The edge in the graph which is the
# substitute is the one with the same terminated end. This can be
# checked rather simply.
#
# Find the edge within k which is the substitute. Because k is a
# 1-arborescence, we know that they is only one such edges
# leading into every vertex.
if len(k.in_edges(e_v, data=weight)) > 1:
raise Exception
sub_u, sub_v, sub_w = next(k.in_edges(e_v, data=weight).__iter__())
k.add_edge(e_u, e_v, **{weight: e_w})
k.remove_edge(sub_u, sub_v)
if (
max(d for n, d in k.in_degree()) <= 1
and len(G) == k.number_of_edges()
and nx.is_weakly_connected(k)
):
# Ascent method calculation
if d[sub_u] == d[e_u] or sub_w == e_w:
# Revert to the original graph
k.remove_edge(e_u, e_v)
k.add_edge(sub_u, sub_v, **{weight: sub_w})
continue
epsilon = (sub_w - e_w) / (d[e_u] - d[sub_u])
if 0 < epsilon < min_epsilon:
min_epsilon = epsilon
# Revert to the original graph
k.remove_edge(e_u, e_v)
k.add_edge(sub_u, sub_v, **{weight: sub_w})
return min_epsilon
# I have to know that the elements in pi correspond to the correct elements
# in the direction of ascent, even if the node labels are not integers.
# Thus, I will use dictionaries to made that mapping.
pi_dict = {}
for n in G:
pi_dict[n] = 0
del n
original_edge_weights = {}
for u, v, d in G.edges(data=True):
original_edge_weights[(u, v)] = d[weight]
dir_ascent, k_d = direction_of_ascent()
while dir_ascent is not None:
max_distance = find_epsilon(k_d, dir_ascent)
for n, v in dir_ascent.items():
pi_dict[n] += max_distance * v
for u, v, d in G.edges(data=True):
d[weight] = original_edge_weights[(u, v)] + pi_dict[u]
dir_ascent, k_d = direction_of_ascent()
nx._clear_cache(G)
# k_d is no longer an individual 1-arborescence but rather a set of
# minimal 1-arborescences at the maximum point of the polytope and should
# be reflected as such
k_max = k_d
# Search for a cycle within k_max. If a cycle exists, return it as the
# solution
for k in k_max:
if len([n for n in k if k.degree(n) == 2]) == G.order():
# Tour found
# TODO: this branch does not restore original_edge_weights of G!
return k.size(weight), k
# Write the original edge weights back to G and every member of k_max at
# the maximum point. Also average the number of times that edge appears in
# the set of minimal 1-arborescences.
x_star = {}
size_k_max = len(k_max)
for u, v, d in G.edges(data=True):
edge_count = 0
d[weight] = original_edge_weights[(u, v)]
for k in k_max:
if (u, v) in k.edges():
edge_count += 1
k[u][v][weight] = original_edge_weights[(u, v)]
x_star[(u, v)] = edge_count / size_k_max
# Now symmetrize the edges in x_star and scale them according to (5) in
# reference [1]
z_star = {}
scale_factor = (G.order() - 1) / G.order()
for u, v in x_star:
frequency = x_star[(u, v)] + x_star[(v, u)]
if frequency > 0:
z_star[(u, v)] = scale_factor * frequency
del x_star
# Return the optimal weight and the z dict
return next(k_max.__iter__()).size(weight), z_star
@nx._dispatchable
def spanning_tree_distribution(G, z):
"""找到Asadpour指数分布的生成树。
使用第7节中的方法解决Asadpour算法[1]_中的最大熵凸规划问题，构建无向生成树的指数分布。
该算法确保生成树中任何边的概率与包含该边的树的概率之和与图中所有生成树的概率之和成比例。
Parameters
----------
G : nx.MultiGraph
用于Held Karp松弛的无向支撑图
z : dict
`held_karp_ascent()` 的输出，Held-Karp解的缩放版本。
Returns
-------
gamma : dict
近似保持 `z` 的边际概率的概率分布。
"""
from math import exp
from math import log as ln
def q(e):
"""在Asadpour论文中，q(e)的值被描述为“边e被包含在一个以概率正比于exp(gamma(T))选择的生成树T中的概率”，这基本上意味着它是该边在整个分布中出现的总概率。
Parameters
----------
e : 元组
描述我们感兴趣的边的 `(u, v)` 元组
Returns
-------
浮点数
根据当前gamma值选择的生成树包含边 `e` 的概率
"""
# Create the laplacian matrices
for u, v, d in G.edges(data=True):
d[lambda_key] = exp(gamma[(u, v)])
G_Kirchhoff = nx.total_spanning_tree_weight(G, lambda_key)
G_e = nx.contracted_edge(G, e, self_loops=False)
G_e_Kirchhoff = nx.total_spanning_tree_weight(G_e, lambda_key)
# Multiply by the weight of the contracted edge since it is not included
# in the total weight of the contracted graph.
return exp(gamma[(e[0], e[1])]) * G_e_Kirchhoff / G_Kirchhoff
# initialize gamma to the zero dict
gamma = {}
for u, v, _ in G.edges:
gamma[(u, v)] = 0
# set epsilon
EPSILON = 0.2
# pick an edge attribute name that is unlikely to be in the graph
lambda_key = "spanning_tree_distribution's secret attribute name for lambda"
while True:
# We need to know that know that no values of q_e are greater than
# (1 + epsilon) * z_e, however changing one gamma value can increase the
# value of a different q_e, so we have to complete the for loop without
# changing anything for the condition to be meet
in_range_count = 0
# Search for an edge with q_e > (1 + epsilon) * z_e
for u, v in gamma:
e = (u, v)
q_e = q(e)
z_e = z[e]
if q_e > (1 + EPSILON) * z_e:
delta = ln(
(q_e * (1 - (1 + EPSILON / 2) * z_e))
/ ((1 - q_e) * (1 + EPSILON / 2) * z_e)
)
gamma[e] -= delta
# Check that delta had the desired effect
new_q_e = q(e)
desired_q_e = (1 + EPSILON / 2) * z_e
if round(new_q_e, 8) != round(desired_q_e, 8):
raise nx.NetworkXError(
f"Unable to modify probability for edge ({u}, {v})"
)
else:
in_range_count += 1
# Check if the for loop terminated without changing any gamma
if in_range_count == len(gamma):
break
# Remove the new edge attributes
for _, _, d in G.edges(data=True):
if lambda_key in d:
del d[lambda_key]
return gamma
[docs]
@nx._dispatchable(edge_attrs="weight")
def greedy_tsp(G, weight="weight", source=None):
"""返回从 `source` 开始的低成本循环及其成本。
这近似解决了旅行商问题。它找到一个包含所有节点的循环，使得销售员可以按顺序访问多个节点，同时最小化总距离。它使用一个简单的贪心算法。本质上，该函数返回一个给定源点的大循环，使得循环的总成本最小化。
Parameters
----------
G : Graph
该图应为一个完整的加权无向图。
所有节点对之间的距离应包含在内。
weight : string, 可选 (默认="weight")
对应边权重的边数据键。
如果任何边没有此属性，则权重设置为1。
source : node, 可选 (默认: 列表中的第一个节点)
起始节点。如果为None，默认为 ``next(iter(G))``
Returns
-------
cycle : list of nodes
返回销售员可以遵循以最小化旅行总权重的循环（节点列表）。
Raises
------
NetworkXError
如果 `G` 不完整，算法会引发异常。
Examples
--------
>>> from networkx.algorithms import approximation as approx
>>> G = nx.DiGraph()
>>> G.add_weighted_edges_from(
... {
... ("A", "B", 3),
... ("A", "C", 17),
... ("A", "D", 14),
... ("B", "A", 3),
... ("B", "C", 12),
... ("B", "D", 16),
... ("C", "A", 13),
... ("C", "B", 12),
... ("C", "D", 4),
... ("D", "A", 14),
... ("D", "B", 15),
... ("D", "C", 2),
... }
... )
>>> cycle = approx.greedy_tsp(G, source="D")
>>> cost = sum(G[n][nbr]["weight"] for n, nbr in nx.utils.pairwise(cycle))
>>> cycle
['D', 'C', 'B', 'A', 'D']
>>> cost
31
Notes
-----
该贪心算法的实现基于以下内容：
- 算法在每次迭代中向解决方案添加一个节点。
- 算法选择一个尚未在循环中的节点，其与前一个节点的连接为循环增加了最小的成本。
贪心算法并不总是给出最佳解决方案。然而，它可以构造一个初始可行解，该解可以作为参数传递给模拟退火或阈值接受等迭代改进算法。
时间复杂度：运行时间为$O(|V|^2)$
"""
# Check that G is a complete graph
N = len(G) - 1
# This check ignores selfloops which is what we want here.
if any(len(nbrdict) - (n in nbrdict) != N for n, nbrdict in G.adj.items()):
raise nx.NetworkXError("G must be a complete graph.")
if source is None:
source = nx.utils.arbitrary_element(G)
if G.number_of_nodes() == 2:
neighbor = next(G.neighbors(source))
return [source, neighbor, source]
nodeset = set(G)
nodeset.remove(source)
cycle = [source]
next_node = source
while nodeset:
nbrdict = G[next_node]
next_node = min(nodeset, key=lambda n: nbrdict[n].get(weight, 1))
cycle.append(next_node)
nodeset.remove(next_node)
cycle.append(cycle[0])
return cycle
[docs]
@py_random_state(9)
@nx._dispatchable(edge_attrs="weight")
def simulated_annealing_tsp(
G,
init_cycle,
weight="weight",
source=None,
temp=100,
move="1-1",
max_iterations=10,
N_inner=100,
alpha=0.01,
seed=None,
):
"""返回旅行商问题的近似解。
该函数使用模拟退火算法来近似通过所有节点的最小成本循环。从次优解开始，模拟退火算法会扰动该解，偶尔接受使解变得更差的改变，以逃离局部最优解。接受这类改变的概率随着迭代次数的增加而降低，以鼓励得到最优结果。简而言之，该函数返回一个从 `source` 开始的循环，使得总成本最小化，并返回该成本。
接受提议改变的概率与一个称为温度（temperature）的参数相关（退火过程在物理上有类似于钢铁冷却硬化的类比）。随着温度的降低，增加成本的移动的概率会下降。
Parameters
----------
G : Graph
`G` 应该是一个完全加权图。
所有节点对之间的距离应该被包含。
init_cycle : 所有节点的列表或 "greedy"
初始解（一个通过所有节点并返回起点的循环）。
此参数没有默认值，以确保你仔细考虑它。
如果为 "greedy"，则使用 `greedy_tsp(G, weight)` 。
其他常见的起始循环是 `list(G) + [next(iter(G))]` 或执行 `threshold_accepting_tsp` 时的 `simulated_annealing_tsp` 的最终结果。
weight : 字符串, 可选 (默认="weight")
对应于边权重的边数据键。
如果任何边没有此属性，则权重设置为1。
source : 节点, 可选 (默认: 列表(G)中的第一个节点)
起始节点。如果为 None，则默认为 ``next(iter(G))``
temp : 整数, 可选 (默认=100)
算法的温度参数。表示初始温度值。
move : "1-1" 或 "1-0" 或 函数, 可选 (默认="1-1")
指示在寻找新试验解时使用何种移动的标志。
字符串表示两种特殊的内置移动：
- "1-1": 1-1 交换，即交换当前解中两个元素的位置。
调用的函数是 :func:`swap_two_nodes` 。
例如，如果在解 ``A = [3, 2, 1, 4, 3]`` 上应用 1-1 交换，
可以通过交换第1和第4个元素得到：
``A' = [3, 2, 4, 1, 3]``
- "1-0": 1-0 交换，即将解中的一个节点移动到一个新位置。
调用的函数是 :func:`move_one_node` 。
例如，如果在解 ``A = [3, 2, 1, 4, 3]`` 上应用 1-0 交换，
可以将第4个元素移动到第2个位置：
``A' = [3, 4, 2, 1, 3]``
你可以提供自己的函数来从一个解移动到邻近解。该函数必须以解和用于控制随机数生成的 `seed` 输入作为参数（参见这里的 `seed` 输入）。你的函数应保持解为一个循环，首尾节点相同且所有其他节点只出现一次。你的函数应返回新解。
max_iterations : 整数, 可选 (默认=10)
当外循环连续迭代次数达到此数值且最佳成本解没有任何变化时，算法宣告完成。
N_inner : 整数, 可选 (默认=100)
内循环的迭代次数。
alpha : 介于 (0, 1) 之间的浮点数, 可选 (默认=0.01)
外循环每次迭代中温度下降的百分比。
seed : 整数, random_state, 或 None (默认)
随机数生成状态的指示符。
参见 :ref:`Randomness<randomness>` 。
Returns
-------
cycle : 节点列表
返回一个旅行商可以遵循的循环（节点列表），以最小化旅行的总权重。
Raises
------
NetworkXError
如果 `G` 不是完全图，算法会引发异常。
Examples
--------
>>> from networkx.algorithms import approximation as approx
>>> G = nx.DiGraph()
>>> G.add_weighted_edges_from(
... {
... ("A", "B", 3),
... ("A", "C", 17),
... ("A", "D", 14),
... ("B", "A", 3),
... ("B", "C", 12),
... ("B", "D", 16),
... ("C", "A", 13),
... ("C", "B", 12),
... ("C", "D", 4),
... ("D", "A", 14),
... ("D", "B", 15),
... ("D", "C", 2),
... }
... )
>>> cycle = approx.simulated_annealing_tsp(G, "greedy", source="D")
>>> cost = sum(G[n][nbr]["weight"] for n, nbr in nx.utils.pairwise(cycle))
>>> cycle
['D', 'C', 'B', 'A', 'D']
>>> cost
31
>>> incycle = ["D", "B", "A", "C", "D"]
>>> cycle = approx.simulated_annealing_tsp(G, incycle, source="D")
>>> cost = sum(G[n][nbr]["weight"] for n, nbr in nx.utils.pairwise(cycle))
>>> cycle
['D', 'C', 'B', 'A', 'D']
>>> cost
31
Notes
-----
模拟退火是一种元启发式局部搜索算法。该算法的主要特点是它接受甚至导致成本增加的解，以逃离低质量的局部最优解。
该算法需要一个初始解。如果没有提供，则由一个简单的贪心算法构造。在每次迭代中，算法会仔细选择一个邻近解。考虑当前解的成本 $c(x)$ 和邻近解的成本 $c(x')$。如果 $c(x') - c(x) <= 0$，则邻近解成为下一次迭代的当前解。否则，算法以概率 $p = exp - ([c(x') - c(x)] / temp)$ 接受邻近解。否则，当前解保持不变。
`temp` 是算法的参数，代表温度。
时间复杂度：
对于内循环的 $N_i$ 次迭代和外循环的 $N_o$ 次迭代，该算法的运行时间为 $O(N_i * N_o * |V|)$。
更多信息和算法的灵感来源请参见：
http://en.wikipedia.org/wiki/Simulated_annealing
"""
if move == "1-1":
move = swap_two_nodes
elif move == "1-0":
move = move_one_node
if init_cycle == "greedy":
# Construct an initial solution using a greedy algorithm.
cycle = greedy_tsp(G, weight=weight, source=source)
if G.number_of_nodes() == 2:
return cycle
else:
cycle = list(init_cycle)
if source is None:
source = cycle[0]
elif source != cycle[0]:
raise nx.NetworkXError("source must be first node in init_cycle")
if cycle[0] != cycle[-1]:
raise nx.NetworkXError("init_cycle must be a cycle. (return to start)")
if len(cycle) - 1 != len(G) or len(set(G.nbunch_iter(cycle))) != len(G):
raise nx.NetworkXError("init_cycle should be a cycle over all nodes in G.")
# Check that G is a complete graph
N = len(G) - 1
# This check ignores selfloops which is what we want here.
if any(len(nbrdict) - (n in nbrdict) != N for n, nbrdict in G.adj.items()):
raise nx.NetworkXError("G must be a complete graph.")
if G.number_of_nodes() == 2:
neighbor = next(G.neighbors(source))
return [source, neighbor, source]
# Find the cost of initial solution
cost = sum(G[u][v].get(weight, 1) for u, v in pairwise(cycle))
count = 0
best_cycle = cycle.copy()
best_cost = cost
while count <= max_iterations and temp > 0:
count += 1
for i in range(N_inner):
adj_sol = move(cycle, seed)
adj_cost = sum(G[u][v].get(weight, 1) for u, v in pairwise(adj_sol))
delta = adj_cost - cost
if delta <= 0:
# Set current solution the adjacent solution.
cycle = adj_sol
cost = adj_cost
if cost < best_cost:
count = 0
best_cycle = cycle.copy()
best_cost = cost
else:
# Accept even a worse solution with probability p.
p = math.exp(-delta / temp)
if p >= seed.random():
cycle = adj_sol
cost = adj_cost
temp -= temp * alpha
return best_cycle
[docs]
@py_random_state(9)
@nx._dispatchable(edge_attrs="weight")
def threshold_accepting_tsp(
G,
init_cycle,
weight="weight",
source=None,
threshold=1,
move="1-1",
max_iterations=10,
N_inner=100,
alpha=0.1,
seed=None,
):
"""返回旅行商问题的近似解。
该函数使用阈值接受方法来近似通过节点的最小成本循环。从次优解开始，阈值接受方法扰动该解，接受任何使解不比增加阈值量更差的改变。接受成本的改进，但也接受导致成本小幅增加的改变。这允许解离开解空间中的次优局部最小值。随着迭代进行，阈值缓慢降低，有助于确保最优解。总之，该函数返回一个从 `source` 开始的循环，其总成本最小化。
Parameters
----------
G : 图
`G` 应该是一个完全加权图。
所有节点对之间的距离应包含在内。
init_cycle : 列表或 "greedy"
初始解（通过所有节点返回起点的循环）。
此参数没有默认值，以确保您考虑它。
如果为 "greedy"，则使用 `greedy_tsp(G, weight)` 。
其他常见的起始循环是 `list(G) + [next(iter(G))]` 或执行 `threshold_accepting_tsp` 时的 `simulated_annealing_tsp` 的最终结果。
weight : 字符串, 可选 (默认="weight")
对应于边权重的边数据键。
如果任何边没有此属性，则权重设置为1。
source : 节点, 可选 (默认: 列表(G)中的第一个节点)
起始节点。如果为 None，则默认为 ``next(iter(G))``
threshold : 整数, 可选 (默认=1)
算法的阈值参数。表示初始阈值的值
move : "1-1" 或 "1-0" 或 函数, 可选 (默认="1-1")
指示在寻找新试验解时使用哪种移动的标志。
字符串表示两种特殊的内置移动：
- "1-1": 1-1 交换，交换当前解中两个元素的位置。
调用的函数是 :func:`swap_two_nodes` 。
例如，如果在解 ``A = [3, 2, 1, 4, 3]`` 中应用 1-1 交换，
我们可以通过交换 1 和 4 元素得到：
``A' = [3, 2, 4, 1, 3]``
- "1-0": 1-0 交换，将解中的一个节点移动到新位置。
调用的函数是 :func:`move_one_node` 。
例如，如果在解 ``A = [3, 2, 1, 4, 3]`` 中应用 1-0 交换，
我们可以将第四个元素转移到第二个位置：
``A' = [3, 4, 2, 1, 3]``
您可以提供自己的函数来从一个解移动到邻近解。该函数必须以解作为输入，并带有控制随机数生成的 `seed` 输入（参见这里的 `seed` 输入）。您的函数应保持解为一个循环，首尾节点相同，其他节点出现一次。您的函数应返回新解。
max_iterations : 整数, 可选 (默认=10)
当外循环连续迭代次数达到此数值且最佳成本解没有任何变化时，声明完成。
N_inner : 整数, 可选 (默认=100)
内循环的迭代次数。
alpha : 介于 (0, 1) 之间的浮点数, 可选 (默认=0.1)
当至少接受一个邻近解时，阈值减少的百分比。
如果没有内循环移动被接受，阈值保持不变。
seed : 整数, random_state, 或 None (默认)
随机数生成状态的指示符。
参见 :ref:`Randomness<randomness>` 。
Returns
-------
cycle : 节点列表
返回旅行商可以遵循的循环（节点列表），以最小化行程的总权重。
Raises
------
NetworkXError
如果 `G` 不是完全图，算法会引发异常。
Examples
--------
>>> from networkx.algorithms import approximation as approx
>>> G = nx.DiGraph()
>>> G.add_weighted_edges_from(
... {
... ("A", "B", 3),
... ("A", "C", 17),
... ("A", "D", 14),
... ("B", "A", 3),
... ("B", "C", 12),
... ("B", "D", 16),
... ("C", "A", 13),
... ("C", "B", 12),
... ("C", "D", 4),
... ("D", "A", 14),
... ("D", "B", 15),
... ("D", "C", 2),
... }
... )
>>> cycle = approx.threshold_accepting_tsp(G, "greedy", source="D")
>>> cost = sum(G[n][nbr]["weight"] for n, nbr in nx.utils.pairwise(cycle))
>>> cycle
['D', 'C', 'B', 'A', 'D']
>>> cost
31
>>> incycle = ["D", "B", "A", "C", "D"]
>>> cycle = approx.threshold_accepting_tsp(G, incycle, source="D")
>>> cost = sum(G[n][nbr]["weight"] for n, nbr in nx.utils.pairwise(cycle))
>>> cycle
['D', 'C', 'B', 'A', 'D']
>>> cost
31
Notes
-----
阈值接受是一种元启发式局部搜索算法。该算法的主要特点是它接受甚至导致成本增加的解，以逃离低质量的局部最优解。
该算法需要一个初始解。该解可以通过一个简单的贪心算法构造。在每次迭代中，它会仔细选择一个邻近解。
考虑当前解的成本 $c(x)$ 和邻近解的成本 $c(x')$。
如果 $c(x') - c(x) <= threshold$，则邻近解成为下一次迭代的当前解，其中阈值称为阈值。
与模拟退火算法相比，阈值接受算法不接受非常低质量的解（由于存在阈值值）。在模拟退火的情况下，即使是非常低质量的解也可以以概率 $p$ 被接受。
时间复杂度：
它的运行时间为 $O(m * n * |V|)$，其中 $m$ 和 $n$ 分别是外循环和内循环运行的次数。
有关算法的更多信息和灵感来源，请参见：
https://doi.org/10.1016/0021-9991(90)90201-B
See Also
--------
simulated_annealing_tsp
"""
if move == "1-1":
move = swap_two_nodes
elif move == "1-0":
move = move_one_node
if init_cycle == "greedy":
# Construct an initial solution using a greedy algorithm.
cycle = greedy_tsp(G, weight=weight, source=source)
if G.number_of_nodes() == 2:
return cycle
else:
cycle = list(init_cycle)
if source is None:
source = cycle[0]
elif source != cycle[0]:
raise nx.NetworkXError("source must be first node in init_cycle")
if cycle[0] != cycle[-1]:
raise nx.NetworkXError("init_cycle must be a cycle. (return to start)")
if len(cycle) - 1 != len(G) or len(set(G.nbunch_iter(cycle))) != len(G):
raise nx.NetworkXError("init_cycle is not all and only nodes.")
# Check that G is a complete graph
N = len(G) - 1
# This check ignores selfloops which is what we want here.
if any(len(nbrdict) - (n in nbrdict) != N for n, nbrdict in G.adj.items()):
raise nx.NetworkXError("G must be a complete graph.")
if G.number_of_nodes() == 2:
neighbor = list(G.neighbors(source))[0]
return [source, neighbor, source]
# Find the cost of initial solution
cost = sum(G[u][v].get(weight, 1) for u, v in pairwise(cycle))
count = 0
best_cycle = cycle.copy()
best_cost = cost
while count <= max_iterations:
count += 1
accepted = False
for i in range(N_inner):
adj_sol = move(cycle, seed)
adj_cost = sum(G[u][v].get(weight, 1) for u, v in pairwise(adj_sol))
delta = adj_cost - cost
if delta <= threshold:
accepted = True
# Set current solution the adjacent solution.
cycle = adj_sol
cost = adj_cost
if cost < best_cost:
count = 0
best_cycle = cycle.copy()
best_cost = cost
if accepted:
threshold -= threshold * alpha
return best_cycle
