中国邮递员问题

中国邮递员问题

个人平常不写算法的博客,主要是太懒太菜,前一段时间离散数学小组展示展示中国邮递员问题现场举例的时候翻车了,这次记录一下中国邮递员问题,来让自己铭记准备不充分理解不透彻队友两行泪

ps:对于我的一些理解错误的地方还请就当看个乐子,还有这次博客巨长(md文件竟然有500+行,对于我平常只写200-300左右,我惊了),测试例子的不多,反正课本上的例子过了,不确定是不是 100% 正确太菜了

中国邮递员问题

问题背景

一个邮递员从邮局出发投递信件,他必须在他所管辖范围内的所有街道至少走一次,最后回到邮局,他自然希望选择一条最短的路线完成投递任务,那么如何选择这样的路线呢?这个问题是中国数学家管梅谷先生首先提出的,因而被称作中国邮递员问题,也可以简称为邮递员问题。

数学建模

要解邮递员问题,首先应将该问题用图来描述。构造无向带权图 G={V,E,W},E 为街道集合,V 中元素为街道的交叉点。街道的长度为该街道对应的边的权,显然所有权均大于0。邮递员问题就变成了求G中一条经过每条边至少一次的回路,使该回路所带权最小的问题,并且称满足以上条件的回路是最优投递路线或最优回路。

理论求解

  • 若G是欧拉图,则最优投递路线为 G 中的任意一条欧拉回路。
  • 若 G 不是欧拉图,则 G 必有奇度顶点,设 G 是连通图,为了消去奇度顶点,必须加若干条重复边,设所得图为 G ,于是求 G 的最优投递路线就等价于求 G 的一条欧拉回路。

加入重复的边需要满足一定的条件:
1) G 的每条边在 G 中至多重复出现一次(否则删掉两条平行边,所有顶点依旧是偶数度,得到了一个总权值更小的欧拉图)
2) G 的每个圈上在 G
中重复出现的边的权之和不超过该圈权的一半。(用通俗点的话来讲就是,一个圆环周长是4,其上两点之间的劣弧长最多是2)

我们可以证明:
设带正权无向连通图 G={V,E,W} ,V’ 为 G 中奇度顶点集,设 |V’|=2k(k≥0),F={e|e∈E 且在求 G 的最优回路时加了重复边},则 F 的导出子图 G[F] 可以表示为以 V’ 中顶点为起点与终点的 k 条不交的最短路径之并。即我们求非欧拉图的中国邮递员问题转换为求 V’ 中顶点为起点和终点的不相交的最短路径的长度

整理一下:
n 为 C 中顶点数。 算法步骤如下:

  1. 若 G 中无奇度顶点,令 G*=G,转2,否则转3。
  2. 求 G* 中的欧拉回路,结束。 -> Fleury 算法
  3. 求 G 中所有奇度顶点对之间的最短路径。 -> Floyd & Dijkstra 算法
  4. 求出奇度顶点之间组合的最小值(最优匹配)
  5. 将 M 中边对应的各最短路径中的边均在 G 中加重复边,得欧拉图G*,转2。

简化问题

问题目前集中到了如何求最小的”组合”,即求奇度顶点集合中顶点为起点和终点的不相交的最短路径的长度之和。顺便也要求出来这个路径经过了那些点,不然就只完成了求最短路径,我们的邮递员还是不知道该怎么走。

算法

暴力

简单易懂,无脑好实现

基本思路是对于奇度顶点,尝试枚举所有的情况的组合来比较得出最小的

例子.png

我们举例来说明:v2,v4,v6,v8 两两组合,可以有限的枚举出组合
{(v2,v4),(v6,v8)},{(v2,v6),(v4,v8)},{(v2,v8),(v4,v6)},我们比较容易比较出来究竟哪个最小。

看起来很舒服,只需要亿个 for 循环就可以搞得定,但如果奇度点增多,导致组合会更加的多,这时如果盲目的枚举会导致算法的时间复杂度垂直上升,这时候就需要一个巧妙的枚举方法:dp 压缩法

经过在网上的学习,我找到了这样的一个算法,它将所有的奇度顶点按照1234重新编号,将它们对应2进制的1234位(从右到左),这样我们可以用1111来表示一个含有这4个点的集合,1101表示含有重新编号过后的含有134这3个点的集合。那么我们对每个集合求最小的不重复的路径的并,最终1111(十进制的15)这个集合的求出的数就是我们想要的。

那应该怎么求呢,我们以上面的例子为例,我们先对 {v2,v4,v6,v8} 重新编号{1,2,3,4}为首先我们定义一个16大小的数组 dp,全部初始化为最大值,对于每次计算一个位置的值的时候,将索引的二进制表示出来,如要计算dp[7](0111),我们就去寻找0111代表的集合{123}的不同组合方式 {12,3}, {13,2}, {23,1},这时我们就可以用 d[3] + d[4] 来更新它了。当然求这个组合是没有任何意义的,这只是让我们能更好的理解这个算法的关键思想:用局部的小值更新大的值。

理解了这个算法的本质时,我们在编写的时候实际并不需要对于dp数组的每个元素都求出其值,我们只需要自己构造就行了,简单来说就虽然也是遍历整个dp数组,但我们每次选择x,y(不在当前的组合中)去构造一个更大的,然后用更大的和我们构造的来比较取小值就可以了,比较的方法就是 min{dp[i]+d[x][y],dp[t]},我们举dp[15]的例子来看,当我们一个2个元素的集合扩充x,y后,比如原本是集合里是{12},加上3和4(对应v6,v8),这时我们比较,如果dp[3]+d[6][8]<dp[15],则更新dp[15]

大概就是这样一个思路,对于初始是3元的集合由于找不到y就完全不需要进行计算,这样相当于我们实际有好多的空间时没有用到的,代码实现如下:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
# -*- coding:utf-8 -*-
import sys
# 方法 1 使用动态规划来搜索最小组合

def floyd(graph_matrix):
n = len(graph_matrix)
# 构造 path 数组
path = []
for i in range(n):
path.append([])
# 初始化
for i in range(n):
for j in range(n):
path[i].append(j)

for k in range(n):
for i in range(n):
for j in range(n):
if graph_matrix[i][k] == sys.maxsize or graph_matrix[k][j] == sys.maxsize:
continue
temp = graph_matrix[i][k] + graph_matrix[k][j]
if temp < graph_matrix[i][j]:
graph_matrix[i][j] = temp
# 更新路径
path[i][j] = path[i][k]
return graph_matrix, path


def shortest_path(graph_matrix, d):
"""
:param graph_matrix: 邻接矩阵
:param d: 度数组
:return: 最短路径
"""
# 度为奇数的顶点
graph_matrix, _ = floyd(graph_matrix)
odd_points = [0] # 初始化包含0便于后续dp数组dp[0]的计算
for i, di in enumerate(d):
if di % 2 == 1:
odd_points.append(i)
# 奇度顶点的个数
odd_num = len(odd_points) - 1
dp = [0] # dp[0] 无一个点所以长度也为0
dp_length = 1 << odd_num
# 初始化数组填入最大的数
for i in range(dp_length):
dp.append(sys.maxsize)
for i in range(dp_length):
x = 1
while (1 << (x-1)) & i:
# 选择不在当前集合的点x
x += 1
for y in range(x+1, odd_num+1): # x+1 <= y <= odd_num
# 选择不在集合内的与x不同的y
if (1 << (y-1)) & i:
continue
else:
# 计算
binary_pt = i | (1 << (x-1)) | (1 << (y-1)) # 新形成的集合
if graph_matrix[odd_points[x]][odd_points[y]] == sys.maxsize:
# x,y 不连通
pass
elif dp[i] == sys.maxsize:
# 对于这个组合无最短路径
pass
else:
trail = dp[i] + graph_matrix[odd_points[x]][odd_points[y]]
# 更新最短路径
dp[binary_pt] = min(dp[binary_pt], trail)
# print(dp[dp_length-1])
return dp[dp_length-1]


def main():
# 记录的所有路求和
n = 0
d = []
min_length = 0
with open("中国邮递员问题测试数据1.txt", 'r') as fp:
scales = fp.readline()
n, m = map(int, scales.split(' '))
# 初始化
for i in range(n):
d.append(0)
graph = [[sys.maxsize]*n for _ in range(n)]
# 读入数据
for line in fp.readlines():
start_point, end_point, weight = map(int, line.split(' '))
# 点标号从 1 开始
start_point -= 1
end_point -= 1
d[start_point] += 1
d[end_point] += 1
# 无向图
graph[start_point][end_point] = weight
graph[end_point][start_point] = weight
min_length += weight
print(graph)
print("最小路径长度为:", end='')
t = 0
for di in d:
if di % 2 == 0:
t += 1
if t == n:
# 所有点读数都是偶数度,为欧拉图
pass
else:
min_length += shortest_path(graph, d)
print(min_length)


if __name__ == '__main__':
main()

KM

方法1虽然好,但是有个问题是求解最小路径是怎么组合出来的相当的让人头疼,不直观,而且毕竟是暴力,奇度顶点多起来的时候数组就开的更加的大

下面介绍的是课本上的标准解法,它通过 km 算法来求出奇度顶点的最小匹配。小问号你是不是有很多的朋友?

我们慢慢来理这个算法,首先我们需要对基础知识进行一些补充的学习:

补基础

二部图

通俗来将就是把一个图分成2部分,这2个部分内部所有点不相连,只有2部分之间有连起来的边,也叫二分图

二部图中的匹配

通俗的描述为从两个部分各取一个点,满足他们之间相连的这样的组合的集合(你可以理解为找男女对象)

完备匹配

通俗的说就是一个匹配 M,M 中包含二部图的所有的点(男女嘉宾全部牵手成功)

最优匹配

将二部图的每条边赋一个权值,对于不连通的部分赋一个0权的边,这就称为了一个完全带权二部图,最优匹配(最佳匹配)就是求一个完备匹配使得边权值总和最大。

增广路径

KM算法用于选择带权二部图的最佳匹配
匈牙利算法用于寻找二分图的最大匹配(完备匹配)
KM 和匈牙利算法都涉及一个很重要的概念:增广路径

增广路径的特点:

  • 有奇数条边
  • 起点在二分图的一侧(x),终点在另一侧(y)
  • 路径的点交替选择x,y侧的点,无重复点
  • 起点终点都是没有匹配的点,其余的点都已经匹配过
  • 路径上所有第奇数边都是目前还没有进入匹配子图的边,第偶数条边都是进入匹配子图中的边

增广路径的性质

  • 加入所有的第奇数条边,删除所有的第偶数条边,匹配数加1

匈牙利算法

匈牙利算法就是一个不断寻找增广路径来使得最终求得的匹配为最大的。

如何找增广路径?
深度搜索,如果xi要和yi匹配,但yi和xk匹配了,就去为xk找新的匹配点(未被匹配过),如果有yk就可以修改匹配方法为(xi,yi),(xk,yk),没有就换一个未匹配的yi与xi尝试匹配。举个例子:

匈牙利算法例子

我们要求一个完备匹配,就对每个点找增广路径,简单来说,为每个点找个匹配,找不到就尝试调整被匹配的点让他换一个,这样2边都能满足了。

比如说14匹配,现在2要和4匹配,就重新为1(和4匹配)找个新的点,5就可以,这样找如果找不到的话就换一个点匹配,假如删去15之间的边,那么2就没法和4匹配,它只能换下一个点(5)尝试匹配

KM 算法

我们回到 KM 这里,KM算法的流程用文字来描述就是:

  1. 对于每个二部图的点,给与一个”顶标”,对于 x 侧,顶标的值记录在 x 数组,第i个元素初始化为 xi 点所有边的最大值,对于y 侧的点,顶标的值初始化为0
  2. 用匈牙利算法在相等子图中求完备匹配(相等子图中的点要满足x[i]+y[j]==w[i][j],即i,j顶标和要为 i,j 之间边的权值)
  3. 如果找不到完备匹配,说明某个点 xi 出发在相等子图中没有找到增广路径。此时需要修改顶标值,对于增广路径上的属于 x 侧的点的顶标减去 d,属于 y 侧的点的顶标加上 d,d 的取值为 min{x[i]+y[j]-w[i][j], i属于增广路,j不属于增广路}
  4. 重复上面直至匈牙利算法求出来完备匹配

可以证明当相等子图的最大匹配为原图的完备匹配时,匹配边的权值和等于所有顶点的顶标和,此匹配即为最佳匹配

而至于为什么要这样修改顶标的值呢?详细的可以看百科,这里直接复制了:

  • 我们一定找到了许多条从 Xi 出发并结束于 X 部的匹配边与未匹配边交替出现的路径,姑且称之为交错树。我们将路径中 X 部的顶点顶标减去一个值 d,交错树中属于 Y 部的顶点顶标加上一个值 d。那么我们会发现:

    • 两端都在交错树中的边 (i,j),其顶标和没有变化。也就是说,它原来属于相等子图,现在仍属于相等子图。
    • 两端都不在交错树中的边 (i,j),其顶标也没有变化。也就是说,它原来属于(或不属于)相等子图,现在仍属于(或不属于)相等子图。
    • X 端不在交错树中,Y 端在交错树中的边 (i,j),它的顶标和会增大。它原来不属于相等子图,现在仍不属于相等子图。
    • X 端在交错树中,Y 端不在交错树中的边 (i,j),它的顶标和会减小。也就说,它原来不属于相等子图,现在可能进入了相等子图,因而使相等子图得到了扩大。
  • 我们修改顶标的目的就是要扩大相等子图。为了保证至少有一条边进入相等子图,我们可以在交错树的边中寻找顶标和与边权之差最小的边,这就是前面说的 d 值。将交错树中属于 X 部的顶点减去 d,交错树中属于 Y 部的顶点加上 d。则可以保证至少有一条边扩充进入相等子图。

举个例子:
KM例子

首先初始化 x 数组和 y 数组(顶标数组),选的都是边的最大值,相等子图是{16, 63, 25} 我们对1找不到完备匹配,得到了一个交错树{163},这时我们计算 d 值,选择不在路径中的 y 点和在路径中的 x 点,可求得d=5,修改顶标后

KM例子

此时相等子图中加入{15}再次寻找,仍无法为3找到匹配,这时我们得到的路径是{152},可求得d=5,再次修改顶标,这次我们会发现{14, 24}加入了,这时我们已经可以找到最大匹配了,为{14,25,36},即为最优匹配

优化

  • 原生的KM为 O(n^4) 的时间复杂度:
    • 需要找O(n)次增广路,每次增广最多需要修改O(n)次顶标,每次修改顶标时由于要枚举边来求 d 值,复杂度为O(n^2)

实际上 KM 算法的复杂度是可以做到O(n^3)的。我们给每个 Y 顶点一个”松弛量” slack,每次开始找增广路时初始化为无穷大。在寻找增广路的过程中,检查边(i,j)时,如果它不在相等子图中,则让slack[j]变成原值与 x[i]+y[j]-w[i][j] 的较小值。这样在修改顶标时,取所有不在交错树中的 Y 顶点的 slack 值中的最小值作为 d 值即可。但还要注意一点:修改顶标后要把所有不在交错树中的 Y 顶点的 slack 值都减去d

代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
# km.py
# 实现的不是太好
class KM:
lx = {} # 一边的标号
ly = {} # 另一边的标号
visit_x = {} # 该点在增广路匹配过程中是否被尝试过
visit_y = {}
match = {} # 匹配的映射关系
slack = {} # slack 优化

def __init__(self, graph: dict):
self.graph = graph
# 初始化均未被访问
for p in graph.keys():
self.visit_x[p] = False
self.lx[p] = -1
self.visit_y[p] = False
self.ly[p] = -1
self.match[p] = -1

def _extended_path(self, x):
# 寻找增广路径
# x 点被选择
self.visit_x[x] = True
for y in self.visit_y.keys():
if self.visit_y[y]:
# 结点被访问过了
continue
if self.graph[x][y] == sys.maxsize:
# 不连通
continue
temp = self.lx[x] + self.ly[y] - self.graph[x][y]
if temp == 0:
# x, y 在相等子图中
self.visit_y[y] = True
if self.match[y] == -1 or self._extended_path(self.match[y]):
# 如果 y 元素没有对应的匹配或者 y 有新的匹配点
self.match[y] = x
return True
elif self.slack[y] > temp:
# slack 优化
self.slack[y] = temp
else:
pass
# 无增广路径
return False

def _solve(self):
for x in self.visit_x:
# 对所有的 x 寻找增广回路
# 重新初始化 slack 中的数值
for y in self.visit_y:
self.slack[y] = sys.maxsize
while True:
# 求新的增广回路,重新初始化访问数组
for y in self.visit_y:
self.visit_y[y] = False
for xt in self.visit_x:
self.visit_x[xt] = False
if self._extended_path(x):
# 有增广回路
break # 跳出循环,判断下一个点
else:
delta = sys.maxsize
# 求得需要对定标进行运算的 d
for y in self.visit_y:
# y 在交错树外,x 在交错树中(dfs失败x一定在交错树中)
if (not self.visit_y[y]) and delta > self.slack[y]:
delta = self.slack[y]
# 对定标的值修正,下一步一定可以再加入一个点
for xt in self.visit_x:
if self.visit_x[xt]:
self.lx[xt] -= delta
for yt in self.visit_y:
if self.visit_y[yt]:
self.ly[yt] += delta
else:
# 修改顶标后,要把所有的slack值都减去delta
# 这是因为x点的标定值减小了delta
# 根据slack的计算也需要变换和x点相关的点的slack值
self.slack[yt] -= delta

def km_solve(self, is_max=True):
# 初始化
if not is_max:
for v_dict in self.graph.values():
for pv in v_dict:
v_dict[pv] = -v_dict[pv]
for x in self.lx:
# 贪心算法, 初始化为点最大的权的边
self.lx[x] = max(list(self.graph[x].values()))
for y in self.ly:
self.ly[y] = 0
self._solve()
part_match = {}
temp_s = 0
for k, v in self.match.items():
temp_s += self.graph[k][v]
if k in part_match.keys() or k in part_match.values():
pass
else:
part_match[k] = v
self.match = part_match
print(part_match)
print(temp_s)
# 匹配整体计算了2次
if is_max:
return temp_s / 2
else:
return -temp_s / 2

一点点细节

到此我们已经可以解决这个问题了,但好像只剩完成了求解长度的部分,那路径呢???

其实这个题目最难的是求回路路径

  • 对于 Floyd 算法,我们求最短的路径时同时要记录节点

  • 适当调整 Fleury 算法的实现:对于带平行边的图的存储,使用邻接矩阵时用a[i][j] = k表示 ij 直接有 k 条平行边,当要删除边时a[i][j]-=1,恢复边的时候a[i][j]+=1

完整代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
# fleury.py
class Fleury:

def __init__(self, graph):
self.graph = graph
self.start = -1
self.route = []
self.num = len(graph) # 图中的点数

def solve(self, current, start):
flag = False # 是否还有与x关联的边
# 倒回会导致上一个点被多记录一次
if self.route and self.route[-1] == current:
pass
else:
self.route.append(current)
for i in range(start, self.num):
# 从 start 开始搜索是否有边
if self.graph[i][current] > 0:
# i 与 current 有边
flag = True
# 删除边
self.graph[i][current] -= 1
self.graph[current][i] -= 1
self.solve(i, 0) # 从 i 开始搜索
break
if not flag:
# 如果没有边与当前节点 current 相连
print(self.route)
self.route.pop()
# if len(self.route) == self.num:
# return self.route
if self.all_zeros():
self.route.append(start)
return self.route
else:
# 回溯
temp = self.route[-1]
# 恢复上一次的 2 条边
self.graph[temp][current] += 1
self.graph[current][temp] += 1
new_start = current + 1
# 路径包含所有的点
self.solve(temp, new_start)

def all_zeros(self):
for i in range(0, self.num):
for j in range(0, self.num):
if self.graph[i][j] != 0:
return False
return True

# km求解中国邮递员.py
import sys
from fleury import Fleury
from km import KM


def euler(index, edge):
f = Fleury(edge)
f.solve(index, index)
return f.route

def floyd(graph_matrix):
n = len(graph_matrix)
# 构造 path 数组
# 如果没有清 0 我们需要手动清0
for i in range(n):
graph_matrix[i][i] = 0
path = []
for i in range(n):
path.append([])
# 初始化
for i in range(n):
for j in range(n):
path[i].append(j)

for k in range(n):
for i in range(n):
for j in range(n):
if graph_matrix[i][k] == sys.maxsize or graph_matrix[k][j] == sys.maxsize:
continue
temp = graph_matrix[i][k] + graph_matrix[k][j]
if temp < graph_matrix[i][j]:
graph_matrix[i][j] = temp
# 更新路径
path[i][j] = path[i][k]

for i in range(n):
graph_matrix[i][i] = sys.maxsize

return graph_matrix, path

d = []
min_length = 0
adjacency_array = [] # 邻接矩阵,用起来完全不是邻接矩阵......不想改了
with open("中国邮递员问题测试数据2.txt", 'r') as fp:
scales = fp.readline()
n, m = map(int, scales.split(' '))
# 初始化
for i in range(n):
d.append(0)
graph = [[sys.maxsize]*n for _ in range(n)]
adjacency_array = [[0]*n for _ in range(n)]
# 读入数据
# 要求读入的数据要满足邻接表的排布
for line in fp.readlines():
start_point, end_point, weight = map(int, line.split(' '))
# 输入的点标号从 1 开始
start_point -= 1
end_point -= 1
# 构造邻接矩阵
adjacency_array[start_point][end_point] += 1
adjacency_array[end_point][start_point] += 1
d[start_point] += 1
d[end_point] += 1
# 无向图
graph[start_point][end_point] = weight
graph[end_point][start_point] = weight
min_length += weight
t = 0
for di in d:
if di % 2 == 0:
t += 1
if t == n:
# 所有点读数都是偶数度,为欧拉图
route = euler(0, adjacency_array)
else:
odd_points = [] # 记录奇数度顶点的数组
for i, di in enumerate(d):
if di % 2 == 1:
odd_points.append(i)
# 奇度顶点的个数
graph_matrix, path = floyd(graph)
# 奇度顶点组成子图
sub_graph = {}
for odd_p in odd_points:
sub_graph[odd_p] = {}
for end_p in odd_points:
sub_graph[odd_p][end_p] = graph_matrix[odd_p][end_p]
km_solution = KM(sub_graph)
# 求最小解
min_length += km_solution.km_solve(False)
# 删去匹配的边
for k, v in km_solution.match.items():
start = k
temp = path[k][v]
while True:
# 存在中转的情况, 以中转为起点
adjacency_array[start][temp] += 1
adjacency_array[temp][start] += 1
if temp == v:
break
else:
pass
start = temp
temp = path[temp][v]

route = euler(0, adjacency_array)
for p in route:
print("({})".format(p+1), end='')
print("最小路径长度为:", end='')
print(min_length)
Author

Ctwo

Posted on

2020-05-04

Updated on

2020-10-25

Licensed under

Comments