Python數(shù)學(xué)建模系列（八）：圖論
1 圖論模型 - Dijkstra
Dijkstra算法能求一個(gè)頂點(diǎn)到另一頂點(diǎn)最短路徑。
基礎(chǔ)概念
無(wú)向圖：若圖中的每條邊都是沒(méi)有方向的，則稱(chēng)該圖為無(wú)向圖。 有向圖：若圖中的每條邊都是有方向的，則稱(chēng)該圖為有向圖。 混合圖：若圖中的部分邊是有方向的，而部分邊是無(wú)方向的，則稱(chēng)該圖為混合圖。
樣例1
如下圖所示，我們需要從①點(diǎn)走到⑨點(diǎn)，每條邊的紅色數(shù) 字代表這條邊的長(zhǎng)度，我們?nèi)绾握业舰俚舰岬淖疃搪窂侥兀?/p>
步驟：
將①標(biāo)記為P，其它標(biāo)記為T(mén)，找出從①出發(fā)當(dāng)前最短的邊所到的點(diǎn)，將該點(diǎn)的T改為P 將所有P點(diǎn)找到可以到達(dá)的T標(biāo)記點(diǎn)上最短的邊，將到達(dá)的點(diǎn)T改為P 重復(fù)步驟，指導(dǎo)終點(diǎn)的T變?yōu)镻
過(guò)程展示：圈加數(shù)字代表每個(gè)頂點(diǎn)，()內(nèi)數(shù)字代表當(dāng)前行走的距離
?1.①(0)
2.①②(1)
3.①②⑤(3)
4.①②⑤(3) ①④(3)
5.①②⑤⑦(5) ①④(3)
6.①②⑤⑦⑥(6) ①④(3)
7.①②⑤⑦⑥(6) ①②⑤③(6) ①④(3)
8.①②⑤⑦⑥(6) ①②⑤③④(7) ①④(3)
9.①②⑤⑦⑥⑧(8) ①②⑤③④(7) ①④(3)
10.①②⑤⑦⑥⑧(8) ①②⑤⑦⑨(8) ①②⑤③④(7) ①④(3)
?
「所以最短的路徑是 ① ② ⑤ ⑦ ⑨， 長(zhǎng)度為8」
帶權(quán)鄰接矩陣：帶權(quán)鄰接矩陣是表示頂點(diǎn)相鄰關(guān)系的矩陣，例如上面那個(gè)圖的帶權(quán)鄰接矩陣如下
?注意：原PPT中有錯(cuò)誤，
?點(diǎn)9 -> 點(diǎn)7 應(yīng)該為inf（上圖已經(jīng)修改）
每個(gè)點(diǎn)和自己的距離為0（主對(duì)角線上元素都是零）
在圖上相鄰兩個(gè)點(diǎn)的如果是連通的，距離就是矩陣的值，無(wú)向 的關(guān)于主對(duì)角線對(duì)稱(chēng)；有向的只有可以過(guò)去的路有數(shù)值
無(wú)法連通的距離就是無(wú)窮，記為inf
Demo代碼
# 運(yùn)行環(huán)境：Vs Code
# PPT中代碼有點(diǎn)錯(cuò)誤 已經(jīng)修改～
# 以下代碼經(jīng)測(cè)試 可行
from collections import defaultdict
from heapq import *
inf = 99999 # 不連通值
mtx_graph = [[0, 1, inf, 3, inf, inf, inf, inf, inf],
[1, 0, 5, inf, 2, inf, inf, inf, inf],
[inf, inf, 0, 1, inf, 6, inf, inf, inf],
[inf, inf, inf, 0, inf, 7, inf, 9, inf],
[inf, 2, 3, inf, 0, 4, 2, inf, 8],
[inf, inf, 6, 7, inf, 0, inf, 2, inf],
[inf, inf, inf, inf, inf, 1, 0, inf, 3],
[inf, inf, inf, inf, inf, inf, 1, 0, 2],
[inf, inf, inf, inf, 8, inf, inf, 2, 0]]
m_n = len(mtx_graph)#帶權(quán)連接矩陣的階數(shù)
edges = [] #保存連通的兩個(gè)點(diǎn)之間的距離(點(diǎn)A、點(diǎn)B、距離)
for i in range(m_n):
for j in range(m_n):
if i!=j and mtx_graph[i][j]!=inf:
edges.append((i,j,mtx_graph[i][j]))
def dijkstra(edges, from_node, to_node):
go_path = []
to_node = to_node-1
g = defaultdict(list)
for l,r,c in edges:
# l:點(diǎn)A r:點(diǎn)B c:距離
# 字典： 點(diǎn)A -> (距離，點(diǎn)B)
g[l].append((c,r))
q, seen = [(0, from_node-1, ())], set()
while q:
(cost, v1, path) = heappop(q)#堆彈出當(dāng)前路徑最小成本
if v1 not in seen:
seen.add(v1)
path = (v1, path)
if v1 == to_node:
break
for c, v2 in g.get(v1, ()):
if v2 not in seen:
heappush(q, (cost+c, v2, path))
if v1!=to_node: #無(wú)法到達(dá)
return float('inf'), []
if len(path)>0:
left=path[0]
go_path.append(left)
right=path[1]
while len(right)>0:
left=right[0]
go_path.append(left)
right=right[1]
go_path.reverse() #逆序變換
for i in range(len(go_path)): #標(biāo)號(hào)加1
go_path[i]=go_path[i]+1
return cost, go_path
leght, path = dijkstra(edges, 1, 9)
print('最短距離為：'+str(leght))
print('前進(jìn)路徑為：'+str(path))
運(yùn)行結(jié)果
最短距離為：8
前進(jìn)路徑為：[1, 2, 5, 7, 9]
2 圖論模型-Floyd
Floyd算法通過(guò)動(dòng)態(tài)規(guī)劃解決任意兩點(diǎn)間的最短路徑（多源最短路徑）的問(wèn)題，「可以正確處理負(fù)權(quán)的最短路徑問(wèn)題」
關(guān)鍵原理：
?枚舉中間點(diǎn)k，找到最小的，作為的最小值
?
關(guān)鍵結(jié)論：
?假設(shè)i和j之間的最短路徑上的結(jié)點(diǎn)集里(不包含i，j)，編號(hào)最大的一個(gè)是x.那么在外循環(huán)k = x時(shí)， 肯定得到了最小值.
?
強(qiáng)歸納法 :
?設(shè)i到x中間編號(hào)最大的是x1，x到j(luò)中間編號(hào)最大的是x2.由于x是i到j(luò)中間編號(hào)最大的，那么顯然
?
根據(jù)結(jié)論，時(shí)候已經(jīng)取得最小值，時(shí)候已經(jīng)取得最 小值
那么時(shí)，和已經(jīng)都取得最小值.
因此時(shí)，執(zhí)行
樣例2
與樣例1相似，只是這次我們「求出從每一個(gè)點(diǎn)到其他點(diǎn)的最短距離和路徑」，每條邊的紅色數(shù)字代表這條邊的長(zhǎng)度。
Demo代碼
from collections import defaultdict
from heapq import *
import numpy as np
inf = 99999 # 不連通值
mtx_graph = [[0, 1, inf, 3, inf, inf, inf, inf, inf],
[1, 0, 5, inf, 2, inf, inf, inf, inf],
[inf, inf, 0, 1, inf, 6, inf, inf, inf],
[inf, inf, inf, 0, inf, 7, inf, 9, inf],
[inf, 2, 3, inf, 0, 4, 2, inf, 8],
[inf, inf, 6, 7, inf, 0, inf, 2, inf],
[inf, inf, inf, inf, inf, 1, 0, inf, 3],
[inf, inf, inf, inf, inf, inf, 1, 0, 2],
[inf, inf, inf, inf, 8, inf, inf, 2, 0]]
def Floyd(graph):
N=len(graph)
A=np.array(graph)
path=np.zeros((N,N))
for i in range(0,N):
for j in range(0,N):
if A[i][j]!=inf:
path[i][j]=j
for k in range(0,N):
for i in range(0,N):
for j in range(0,N):
# 原PPT這一句代碼有錯(cuò) 寫(xiě)成了 A[i][j]+A[k][j]<A[i][j]
# 正確應(yīng)該是：A[i][k]+A[k][j]<A[i][j]
if A[i][k]+A[k][j]<A[i][j]:
A[i][j]=A[i][k]+A[k][j]
path[i][j]=path[i][k]
for i in range(0,N):
for j in range(0,N):
path[i][j]=path[i][j]+1
print('距離 = ')
print(A)
print('路徑 = ')
print(path)
Floyd(mtx_graph)
運(yùn)行結(jié)果
?注意：原PPT 結(jié)果中 最后一行有錯(cuò) 應(yīng)該是
?5. 5. 8. 8. 5. 8. 8. 8. 9.PPT中是：5. 5. 7.7. 5. 7. 7. 8. 9.因?yàn)樵谧畛醯?mtx_graph設(shè)置就錯(cuò)了：點(diǎn)9到點(diǎn)7 應(yīng)該是inf 而不是 3
結(jié)論：
從點(diǎn)1 到 點(diǎn)9 的距離為8（結(jié)果中“距離”矩陣中第一行第九列數(shù)值） 路徑為【1 - 2 - 5 - 7 - 9】
怎么看路徑呢？
從最后結(jié)果中查看”路徑“矩陣 點(diǎn)1 到 點(diǎn)9 首先看第一行第九列 數(shù)值為2 這是中間位置2，再看第二行第九列，數(shù)值為5 再看第五行第九列，數(shù)值為7 再看第七行第九列，數(shù)值為9 得到最終路徑【1 - 2 - 5 - 7 - 9】
3 機(jī)場(chǎng)航線設(shè)計(jì)
數(shù)據(jù)集來(lái)自航空業(yè)，有一些關(guān)于航線的基本信息。有某段旅程的起始點(diǎn)和目的地。還有一些列表示每段旅程的到達(dá)和起飛時(shí)間。這個(gè)數(shù)據(jù)集非常適合作為圖進(jìn)行分析。想象一下通過(guò)航線（邊）連接的幾個(gè)城市（節(jié)點(diǎn)）。
如果你是航空公司，你可以問(wèn)如下幾個(gè)問(wèn)題：
從A到B的最短途徑是什么？分別從距離和時(shí)間角度考慮。 有沒(méi)有辦法從C到D？ 哪些機(jī)場(chǎng)的交通最繁忙？ 哪個(gè)機(jī)場(chǎng)位于大多數(shù)其他機(jī)場(chǎng)“之間”？這樣它就可以變成當(dāng)?shù)氐囊粋€(gè)中轉(zhuǎn)站。
0、Airlines.csv數(shù)據(jù)
1、數(shù)據(jù)導(dǎo)入、觀察變量
import numpy as np
import pandas as pd
data = pd.read_csv('../Profile/Airlines.csv')
data.shape #數(shù)據(jù)大小
# (100, 16) 100行16列
運(yùn)行結(jié)果
「查看各個(gè)變量的類(lèi)型」
data.dtypes # 各變量的類(lèi)型
運(yùn)行結(jié)果
2、數(shù)據(jù)清洗
#將sched_dep_time轉(zhuǎn)換為'std'—預(yù)定的出發(fā)時(shí)間
data['std'] = data.sched_dep_time.astype(str).str.replace('(\d{2}$)', '') + ':' + data.sched_dep_time.astype(str).str.extract('(\d{2}$)', expand=False) + ':00'
#將sched_arr_time轉(zhuǎn)換為“sta”—預(yù)定到達(dá)時(shí)間
data['sta'] = data.sched_arr_time.astype(str).str.replace('(\d{2}$)', '') + ':' + data.sched_arr_time.astype(str).str.extract('(\d{2}$)', expand=False) + ':00'
#將dep_time轉(zhuǎn)換為'atd' -實(shí)際出發(fā)時(shí)間
data['atd'] = data.dep_time.fillna(0).astype(np.int64).astype(str).str.replace('(\d{2}$)', '') + ':' + data.dep_time.fillna(0).astype(np.int64).astype(str).str.extract('(\d{2}$)', expand=False) + ':00'
#將arr_time轉(zhuǎn)換為'ata' -實(shí)際到達(dá)時(shí)間
data['ata'] = data.arr_time.fillna(0).astype(np.int64).astype(str).str.replace('(\d{2}$)', '') + ':' + data.arr_time.fillna(0).astype(np.int64).astype(str).str.extract('(\d{2}$)', expand=False) + ':00'
3、時(shí)間信息合并
data['date'] = pd.to_datetime(data[['year', 'month', 'day']])
data = data.drop(columns = ['year', 'month', 'day'])
4、創(chuàng)建圖
import networkx as nx
FG = nx.from_pandas_edgelist(data, source='origin', target='dest', edge_attr=True,)
# 查看所有節(jié)點(diǎn)
FG.nodes()
運(yùn)行結(jié)果
#NodeView(('EWR', 'MEM', 'LGA', 'FLL', 'SEA', 'JFK', 'DEN', 'ORD', 'MIA', 'PBI', 'MCO', 'CMH', 'MSP', 'IAD', 'CLT', 'TPA', 'DCA', 'SJU', 'ATL', 'BHM', 'SRQ', 'MSY', 'DTW', 'LAX', 'JAX', 'RDU', 'MDW', 'DFW', 'IAH', 'SFO', 'STL', 'CVG', 'IND', 'RSW', 'BOS', 'CLE'))
「查詢圖中的所有邊」
# 查看所有邊
FG.edges()
運(yùn)行結(jié)果
#EdgeView([('EWR', 'MEM'), ('EWR', 'SEA'), ('EWR', 'MIA'), ('EWR', 'ORD'), ('EWR', 'MSP'), ('EWR', 'TPA'), ('EWR', 'MSY'), ('EWR', 'DFW'), ('EWR', 'IAH'), ('EWR', 'SFO'), ('EWR', 'CVG'), ('EWR', 'IND'), ('EWR', 'RDU'), ('EWR', 'IAD'), ('EWR', 'RSW'), ('EWR', 'BOS'), ('EWR', 'PBI'), ('EWR', 'LAX'), ('EWR', 'MCO'), ('EWR', 'SJU'), ('LGA', 'FLL'), ('LGA', 'ORD'), ('LGA', 'PBI'), ('LGA', 'CMH'), ('LGA', 'IAD'), ('LGA', 'CLT'), ('LGA', 'MIA'), ('LGA', 'DCA'), ('LGA', 'BHM'), ('LGA', 'RDU'), ('LGA', 'ATL'), ('LGA', 'TPA'), ('LGA', 'MDW'), ('LGA', 'DEN'), ('LGA', 'MSP'), ('LGA', 'DTW'), ('LGA', 'STL'), ('LGA', 'MCO'), ('LGA', 'CVG'), ('LGA', 'IAH'), ('FLL', 'JFK'), ('SEA', 'JFK'), ('JFK', 'DEN'), ('JFK', 'MCO'), ('JFK', 'TPA'), ('JFK', 'SJU'), ('JFK', 'ATL'), ('JFK', 'SRQ'), ('JFK', 'DCA'), ('JFK', 'DTW'), ('JFK', 'LAX'), ('JFK', 'JAX'), ('JFK', 'CLT'), ('JFK', 'PBI'), ('JFK', 'CLE'), ('JFK', 'IAD'), ('JFK', 'BOS')])
「計(jì)算圖的平均邊密度」
#注意我們的100行數(shù)據(jù)都來(lái)自3個(gè)機(jī)場(chǎng)
nx.draw_networkx(FG, with_labels=True)
nx.algorithms.degree_centrality(FG)
# 圖的平均邊密度
nx.density(FG)
# 0.09047619047619047
運(yùn)行結(jié)果（對(duì)圖進(jìn)行可視化）
「求圖中所有路徑的平均最短路徑長(zhǎng)度」
#圖中所有路徑的平均最短路徑長(zhǎng)度
nx.average_shortest_path_length(FG)
# 2.36984126984127
運(yùn)行結(jié)果
「對(duì)于一個(gè)度為k的節(jié)點(diǎn)-求的鄰居度的平均值」
#對(duì)于一個(gè)度為k的節(jié)點(diǎn)-它的鄰居度的平均值是多少
nx.average_degree_connectivity(FG) # For a node of degree k - What is the average of its neighbours' degree?
#{20: 1.95, 1: 19.307692307692307, 2: 19.0625, 17: 2.0588235294117645, 3: 19.0}
運(yùn)行結(jié)果
5、計(jì)算航線
# 找出所有 JAX 到 DFW 的航線
for path in nx.all_simple_paths(FG, source='JAX', target='DFW'):
print(path)
#找出從JAX到DFW的dijkstra路徑。
dijpath = nx.dijkstra_path(FG, source='JAX', target='DFW')
dijpath # ['JAX', 'JFK', 'SEA', 'EWR', 'DFW'] 下圖中的最后一行
運(yùn)行結(jié)果
「從所有航線中找出 飛行時(shí)間最短的路徑」
# 從所有航線中找出 飛行時(shí)間最短的路徑
shortpath = nx.dijkstra_path(FG, source='JAX', target='DFW', weight='air_time')
shortpath
# ['JAX', 'JFK', 'BOS', 'EWR', 'DFW']
運(yùn)行結(jié)果：
注意事項(xiàng)
起始點(diǎn)和目的地可以作為節(jié)點(diǎn)，其他信息應(yīng)當(dāng)作為節(jié)點(diǎn)或邊屬性；單條邊可以被認(rèn)為是一段旅程。這樣的旅程將有不同的時(shí)間，航班號(hào)，飛機(jī)尾號(hào)等相關(guān)信息。 注意到年，月，日和時(shí)間信息分散在許多列；想創(chuàng)建一 個(gè)包含所有這些信息的日期時(shí)間列，還需要將預(yù)計(jì)的 （scheduled）和實(shí)際的(actual)到達(dá)離開(kāi)時(shí)間分開(kāi)；最終 應(yīng)該有4個(gè)日期時(shí)間列（預(yù)計(jì)到達(dá)時(shí)間、預(yù)計(jì)起飛時(shí)間、 實(shí)際到達(dá)時(shí)間和實(shí)際起飛時(shí)間） 時(shí)間格式問(wèn)題 數(shù)據(jù)類(lèi)型問(wèn)題 NaN值的麻煩
結(jié)語(yǔ)
學(xué)習(xí)來(lái)源：B站及其課堂PPT，對(duì)其中代碼進(jìn)行了復(fù)現(xiàn)
?https://www.bilibili.com/video/BV12h411d7Dm
PPT中錯(cuò)誤有點(diǎn)多 用了挺久的時(shí)間學(xué)習(xí)原理找錯(cuò)誤 哎～
以上代碼均已運(yùn)行成功！環(huán)境：Vs Code
?
「文章僅作為學(xué)習(xí)筆記，記錄從0到1的一個(gè)過(guò)程」
希望對(duì)您有所幫助，如有錯(cuò)誤歡迎小伙伴指正～