import random as rd
import copy
import time
from pyscipopt import Model, quicksum
def limit_generate(nums, n_min, n_max):
limit = [[int((n_max-n_min) * rd.random() + n_min) if i != j else 0 for j in range(nums)] for i in range(nums)]
for i in range(nums-1):
for j in range(i+1, nums):
limit[i][j] = limit[j][i]
return limit
def sub_LL_print(L, digit):
s = "["
for i in range(len(L)):
s = s + str(L[i])
temp = L[i]
digit_temp = 1
while temp//10>0:
digit_temp += 1
temp = temp//10
for j in range(digit-digit_temp):
s += " "
if i < len(L)-1:
s += ", "
s += "]"
return s
def LL_print(name, LL):
LL1 = [[int(LL[i][j]) for j in range(len(LL[i]))] for i in range(len(LL))]
num_max = max([max(LL1[i]) for i in range(len(LL1))])
digit = 1
while num_max//10>0:
digit += 1
num_max = num_max // 10
print("\n"+name+":")
print("["+sub_LL_print(LL[0], digit)+",")
for i in range(1,len(LL)-1):
print(" "+sub_LL_print(LL[i], digit)+",")
print(" " + sub_LL_print(LL[-1], digit) + "]")
def sub_LL_print_str(L, digit):
s = "["
for i in range(len(L)):
s = s + L[i]
temp = L[i]
digit_temp = len(L[i])
for j in range(digit-digit_temp):
s += " "
if i < len(L)-1:
s += ", "
s += "]"
return s
def LL_print_str(name, LL):
digit = max([max([len(LL[i][j]) for j in range(len(LL[i]))]) for i in range(len(LL))])
print("\n",name+":")
print("["+sub_LL_print_str(LL[0], digit)+",")
for i in range(1,len(LL)-1):
print(" "+sub_LL_print_str(LL[i], digit)+",")
print(" " + sub_LL_print_str(LL[-1], digit) + "]")
def LF_generate(F):
n = len(F)
LF = [F]
Fi = F
for i in range(n-2):
Fi = [[1 if sum([Fi[j][l]*F[l][k] for l in range(n)]) > 0 else 0 for k in range(n)] for j in range(n)]
LF.append(copy.deepcopy(Fi))
return LF
def A_matrix(limit):
n = len(limit)
L = copy.deepcopy(limit)
F = [[1 if L[i][j] > 0 else 0 for j in range(n)] for i in range(n)]
LF = LF_generate(F)
R = []
N = []
while sum([LF[i][0][n-1] for i in range(len(LF))]) > 0:
for i in range(len(LF)):
if LF[i][0][n-1] == 0:
continue
#说明是i阶邻居，依次找到i-1,...,1阶邻居
route = [n-1]
nj = n-1
for j in range(i-1,-1,-1):
Fj = LF[j]
for k in range(n):
if Fj[0][k] > 0 and F[k][nj] > 0:
nj = k
route.append(nj)
break
route.append(0)
route1 = [route[len(route)-1-j] for j in range(len(route))]
break
flow = min([L[route1[i]][route1[i+1]] for i in range(len(route1)-1)])
R.append(route1)
N.append(flow)
for i in range(len(route1)-1):
L[route1[i]][route1[i+1]] -= flow
F = [[1 if L[i][j] > 0 else 0 for j in range(n)] for i in range(n)]
LF = LF_generate(F)
F1 = [[limit[i][j] - L[i][j] for j in range(n)] for i in range(n)]
return F1, R, N
def sub_IP(L, time_limit):
n = len(L)
M = max([1 / L[i][j] if L[i][j] > 0 else 0 for i in range(n) for j in range(n)] + [1])
m = min([1 / L[i][j] if L[i][j] > 0 else 1 for i in range(n) for j in range(n)] + [1])
model = Model("Maximum_flow_problem_sub_IP")
F = [[model.addVar(vtype="I", name="F[%s,%s]" % (i, j)) for j in range(n)] for i in range(n)]
I = [[model.addVar(vtype="B", name="I[%s,%s]" % (i, j)) for j in range(n)] for i in range(n)]
model.setObjective(quicksum(F[0][j] for j in range(n)), "maximize")
#变量大于等于0且小于限制
for i in range(n):
for j in range(n):
model.addCons(F[i][j] - L[i][j] <= 0)
model.addCons(F[i][j] >= 0)
#其他点进出平衡
for i in range(1,n-1):
model.addCons(quicksum(F[i][j] for j in range(n)) - quicksum(F[j][i] for j in range(n)) == 0)
#起点不进、终点不出
for i in range(n):
model.addCons(F[i][0] == 0)
model.addCons(F[n-1][i] == 0)
# flag约束
for i in range(n):
for j in range(n):
model.addCons(I[i][j] - M * F[i][j] <= 0)
model.addCons(I[i][j] - m * F[i][j] >= 0)
model.addCons(I[i][j] + I[j][i] <= 1)
model.setRealParam("limits/time", time_limit)
model.optimize()
print("\ngap:", model.getGap())
# 拿结果
F1 = [[round(model.getVal(F[i][j])) for j in range(n)] for i in range(n)]
c = sum(F1[0])
return c
def IP(L, time_limit):
c = sub_IP(L, time_limit)
print("\nc:", c)
n = len(L)
M = max([1/L[i][j] if L[i][j] > 0 else 0 for i in range(n) for j in range(n)] + [1])
m = min([1/L[i][j] if L[i][j] > 0 else 1 for i in range(n) for j in range(n)] + [1])
model = Model("Maximum_flow_problem_IP")
F = [[model.addVar(vtype="I", name="F[%s,%s]" % (i, j)) for j in range(n)] for i in range(n)]
I = [[model.addVar(vtype="B", name="I[%s,%s]" % (i, j)) for j in range(n)] for i in range(n)]
model.setObjective(quicksum(F[i][j] for i in range(n) for j in range(n)), "minimize")
#流量不少
model.addCons(quicksum(F[0][j] for j in range(n)) == c)
#变量大于等于0且小于限制
for i in range(n):
for j in range(n):
model.addCons(F[i][j] - L[i][j] <= 0)
model.addCons(F[i][j] >= 0)
#其他点进出平衡
for i in range(1,n-1):
model.addCons(quicksum(F[i][j] for j in range(n)) - quicksum(F[j][i] for j in range(n)) == 0)
#起点不进、终点不出
for i in range(n):
model.addCons(F[i][0] == 0)
model.addCons(F[n-1][i] == 0)
# flag约束
for i in range(n):
for j in range(n):
model.addCons(I[i][j] - M * F[i][j] <= 0)
model.addCons(I[i][j] - m * F[i][j] >= 0)
model.addCons(I[i][j] + I[j][i] <= 1)
model.setRealParam("limits/time", time_limit)
model.optimize()
print("\ngap:", model.getGap())
# 拿结果
F1 = [[round(model.getVal(F[i][j])) for j in range(n)] for i in range(n)]
return F1
def sub_IP2(L, time_limit):
M = sum([sum(L[i]) for i in range(len(L))])
n = len(L)
A = [[0 for j in range(n*n)] for i in range(n)]
for i in range(n):
for j in range(n):
A[i][i*n+j] += 1
A[j][i*n+j] -= 1
# 流量限制
b = [L[i][j] for i in range(n) for j in range(n)]
# 终点到起点限制调大
b[(n - 1) * n] = 10 * M
model = Model("Maximum_flow_problem_IP2")
F = [model.addVar(vtype="I", name="F[%s,%s]" % (i, j)) for i in range(n) for j in range(n)]
model.setObjective(F[(n-1)*n], "maximize")
#流大于等于0
for i in range(n*n):
model.addCons(F[i] >= 0)
#进出平衡
for i in range(n):
model.addCons(quicksum(A[i][j]*F[j] for j in range(n*n)) == 0)
for i in range(n*n):
model.addCons(F[i] - b[i] <= 0)
model.setRealParam("limits/time", time_limit)
model.optimize()
print("\ngap:", model.getGap())
# 拿结果
c = round(model.getVal(F[(n-1)*n]))
return c
def IP2(L, time_limit):
c = sub_IP2(L, time_limit)
print("\nc:", c)
M = sum([sum(L[i]) for i in range(len(L))])
n = len(L)
A = [[0 for j in range(n*n)] for i in range(n)]
for i in range(n):
for j in range(n):
A[i][i*n+j] += 1
A[j][i*n+j] -= 1
# 流量限制
b = [L[i][j] for i in range(n) for j in range(n)]
# 终点到起点限制调大
b[(n - 1) * n] = 10 * M
model = Model("Maximum_flow_problem_IP2")
F = [model.addVar(vtype="I", name="F[%s,%s]" % (i, j)) for i in range(n) for j in range(n)]
model.setObjective(quicksum(F[i] for i in range(n*n)), "minimize")
# 流量不少
model.addCons(F[(n-1)*n] == c)
#流大于等于0
for i in range(n*n):
model.addCons(F[i] >= 0)
#进出平衡
for i in range(n):
model.addCons(quicksum(A[i][j]*F[j] for j in range(n*n)) == 0)
for i in range(n*n):
model.addCons(F[i] - b[i] <= 0)
model.setRealParam("limits/time", time_limit)
model.optimize()
print("\ngap:", model.getGap())
# 拿结果
F1 = [[round(model.getVal(F[i*n+j])) for j in range(n)] for i in range(n)]
F1[n-1][0] = 0
return F1
def sub_LL_print(L, digit):
s = "["
for i in range(len(L)):
s = s + str(L[i])
temp = L[i]
digit_temp = 1
while temp//10>0:
digit_temp += 1
temp = temp//10
for j in range(digit-digit_temp):
s += " "
if i < len(L)-1:
s += ", "
s += "]"
return s
def Route_print(name, R, N):
num_max = max([max(R[i]) for i in range(len(R))])
digit_max = 1
while num_max // 10 > 0:
digit_max += 1
num_max = num_max // 10
R1 = []
print("\n" + name + ":")
for i in range(len(R)):
r = ""
for j in range(len(R[i])):
num = R[i][j]
digit = 1
while num // 10 > 0:
digit += 1
num = num // 10
r = r + str(R[i][j])
for k in range(digit_max-digit):
r = r + " "
if j < len(R[i])-1:
r = r + "-"
r = r + ": " + str(N[i])
print(r)
def Ford_Fulkerson_opt(limit):
#0-初始化点、边、可行流
M = 10000
n = len(limit)
edges = {}
for i in range(n):
temp = [[i,j,limit[i][j],0] for j in range(n) if i != j and limit[i][j] > 0]
edges[i] = temp
edges1 = {}
for i in range(n):
temp = [[j,i,limit[j][i],0] for j in range(n) if i != j and limit[j][i] > 0]
edges1[i] = temp
flags = [0 for i in range(n)]
#1
flags[0] = [0,1,M]
#已标未查集
L = [0]
#已标已查集
S = []
idex = 0
while idex < len(L):
#2
u = L[idex]
for edge in edges[u]:
v,l,used = edge[1],edge[2],edge[3]
#如果未被标号且可以可增
if flags[v] == 0 and used < l:
flags[v] = [u, 1, min(flags[u][2], l-used)]
L.append(v)
for edge in edges1[u]:
v,l,used = edge[0],edge[2],edge[3]
# 如果未被标号且可以有流
if flags[v] == 0 and used > 0:
flags[v] = [u, -1, used]
L.append(v)
#3
S.append(L[idex])
del L[idex]
if n-1 not in L:
#4
if L == []:
break
else:
continue
else:
#5
z = n-1
#6
flagz = flags[z]
w = flagz[0]
flow = flagz[2]
while z != 0:
if flagz[1] == 1:
for i in range(len(edges[w])):
if z == edges[w][i][1]:
edges[w][i][3] += flow
break
for i in range(len(edges1[z])):
if w == edges1[z][i][0]:
edges1[z][i][3] += flow
break
else:
for i in range(len(edges[z])):
if w == edges[z][i][1]:
edges[z][i][3] -= flow
break
for i in range(len(edges1[w])):
if z == edges1[w][i][0]:
edges1[w][i][3] -= flow
break
z = w
flagz = flags[z]
w = flagz[0]
flags = [0 for i in range(n)]
flags[0] = [0, 1, M]
L = [0]
S = []
idex = 0
F = [[0 for j in range(n)] for i in range(n)]
for i in edges:
for edge in edges[i]:
F[edge[0]][edge[1]] = edge[3]
return F
def NF_generate(F, limit):
n = len(limit)
NF = [[0 for j in range(n)] for i in range(n)]
for i in range(n):
for j in range(n):
if limit[i][j] > 0:
NF[i][j] += limit[i][j] - F[i][j]
NF[j][i] += F[i][j]
return NF
def Vs_generate(NF, F, limit):
#初始化边集
n = len(NF)
edges = {}
for i in range(n):
temp = [[i,j,limit[i][j], F[i][j]] for j in range(n) if i != j and limit[i][j] > 0]
edges[i] = temp
edges1 = {}
for i in range(n):
temp = [[j, i, limit[j][i], F[j][i]] for j in range(n) if i != j and limit[j][i] > 0]
edges1[i] = temp
V_sum = [0]
Vs = [[0]]
#节点分层
while True:
temp = []
for p in Vs[-1]:
for i in range(len(edges[p])):
q = edges[p][i][1]
if q not in V_sum and edges[p][i][3] < edges[p][i][2]:
temp.append(q)
V_sum.append(q)
for i in range(len(edges1[p])):
q = edges1[p][i][0]
if q not in V_sum and edges1[p][i][3] > 0:
temp.append(q)
V_sum.append(q)
if len(temp) == 0:
break
Vs.append(temp)
return Vs
def AN_generate(Vs, NF):
n = len(NF)
AN = copy.deepcopy(NF)
idex_y = -1
for i in range(len(Vs)):
if n-1 in Vs[i]:
idex_y = i
#若汇点不在分层网络中，则以得到最大流
if idex_y == -1:
return -1, [], []
#删去比y层数高的顶点和同层的顶点，删去从高层指向低层的弧及同层顶点拣的弧
for i in range(idex_y, len(Vs)):
for p in Vs[i]:
if p != n-1:
for j in range(n):
AN[p][j] = 0
AN[j][p] = 0
for i in range(idex_y):
for j in range(i+1, idex_y+1):
for u in Vs[i]:
for v in Vs[j]:
AN[v][u] = 0
for i in range(idex_y):
for u in Vs[i]:
for v in Vs[i]:
AN[v][u] = 0
#改进分层网络
idex = len(Vs)-1
while idex > 0:
vs = Vs[idex]
us = copy.deepcopy(Vs[idex-1])
i = 0
while i < len(us):
if sum([AN[us[i]][vs[j]] for j in range(len(vs))]) <= 0:
for j in range(n):
AN[us[i]][j] = 0
AN[j][us[i]] = 0
del us[i]
else:
i += 1
Vs[idex - 1] = us
idex -= 1
return 1, AN, Vs
#一次只加一条增广路
def F_iteration(F, AN):
n = len(F)
edges = {}
for i in range(n):
temp = [[i,j,AN[i][j],0] for j in range(n) if i != j and AN[i][j] > 0]
edges[i] = temp
R = [0]
while True:
p = R[-1]
edges_p = edges[p]
if len(edges_p) == 0:
u = R[-2]
for i in range(len(edges[u])):
if p == edges[u][i][1]:
del edges[u][i]
del R[-1]
break
continue
R.append(edges_p[0][1])
if edges_p[0][1] == n-1:
break
flow = min([AN[R[i]][R[i+1]] for i in range(len(R)-1)])
F1 = copy.deepcopy(F)
for i in range(len(R)-1):
u,v = R[i], R[i+1]
if F1[v][u] >= flow:
F1[v][u] -= flow
else:
F1[u][v] += flow - F1[v][u]
F1[v][u] = 0
return F1
#深度优先迭代，一次可增加多条增广路
def F_iteration_DFS(F, AN):
n = len(F)
edges = {}
for i in range(n):
temp = [[i,j,AN[i][j],0] for j in range(n) if i != j and AN[i][j] > 0]
edges[i] = temp
R = [0]
while True:
p = R[-1]
edges_p = edges[p]
if len(edges_p) == 0:
#起点没有可以出去的边了
if len(R) <= 1:
break
u = R[-2]
for i in range(len(edges[u])):
if p == edges[u][i][1]:
del edges[u][i]
del R[-1]
break
continue
R.append(edges_p[0][1])
if edges_p[0][1] == n-1:
flow = min([AN[R[i]][R[i+1]] for i in range(len(R)-1)])
F1 = copy.deepcopy(F)
for i in range(len(R)-1):
u,v = R[i], R[i+1]
if F1[v][u] >= flow:
F1[v][u] -= flow
else:
F1[u][v] += flow - F1[v][u]
F1[v][u] = 0
if AN[u][v] >= flow:
AN[u][v] -= flow
else:
AN[u][v] = 0
AN[v][u] += flow - AN[u][v]
F = F1
edges = {}
for i in range(n):
temp = [[i, j, AN[i][j], 0] for j in range(n) if i != j and AN[i][j] > 0]
edges[i] = temp
R = [0]
return F
#推拉流
def F_iteration_preflow_push(F, AN, Vs):
n = len(F)
Vs1 = copy.deepcopy(Vs)
for i in range(len(Vs1)):
if n-1 in Vs1[i]:
Vs1[i] = [n-1]
del Vs1[i+1:]
break
#根据点的进出能力确定该点能力
points = []
T = []
for i in range(n):
if sum(AN[i]) > 0 or sum([AN[j][i] for j in range(n)]) > 0:
points.append(i)
if i == 0:
T.append(sum(AN[i]))
elif i == n-1:
T.append(sum([AN[j][i] for j in range(n)]))
else:
T.append(min(sum(AN[i]), sum([AN[j][i] for j in range(n)])))
flag_temp = 0
idex_min = points[0]
for i in range(len(points)):
if T[i] > 0:
if flag_temp == 0:
flag_temp = 1
idex_min = i
elif T[idex_min] > T[i]:
idex_min = i
level_min = 0
for i in range(len(Vs1)):
if points[idex_min] in Vs1[i]:
level_min = i
break
dF = [[0 for j in range(n)] for i in range(n)]
#往后流
nums = [T[idex_min] if i == points[idex_min] else 0 for i in Vs1[level_min]]
idex = level_min
while idex < len(Vs1)-1:
nums1 = [0 for i in range(len(Vs1[idex+1]))]
for i in range(len(Vs1[idex])):
left = nums[i]
u = Vs1[idex][i]
for j in range(len(Vs1[idex+1])):
if left == 0:
break
v = Vs1[idex+1][j]
if dF[u][v] < AN[u][v]:
if dF[u][v] + left <= AN[u][v]:
dF[u][v] += left
nums1[j] += left
left = 0
else:
left -= (AN[u][v] - dF[u][v])
nums1[j] += (AN[u][v] - dF[u][v])
dF[u][v] = AN[u][v]
nums = nums1
idex += 1
#往前推
nums = [T[idex_min] if i == points[idex_min] else 0 for i in Vs1[level_min]]
idex = level_min
while idex > 0:
nums1 = [0 for i in range(len(Vs1[idex-1]))]
for i in range(len(Vs1[idex])):
left = nums[i]
v = Vs1[idex][i]
for j in range(len(Vs1[idex-1])):
if left == 0:
break
u = Vs1[idex-1][j]
if dF[u][v] < AN[u][v]:
if dF[u][v] + left <= AN[u][v]:
dF[u][v] += left
nums1[j] += left
left = 0
else:
left -= (AN[u][v] - dF[u][v])
nums1[j] += (AN[u][v] - dF[u][v])
dF[u][v] = AN[u][v]
nums = nums1
idex -= 1
return [[F[i][j] + dF[i][j] for j in range(n)] for i in range(n)]
def dinic(limit, type_push):
#初始化可行流及增量网络
n = len(limit)
F = [[0 for j in range(n)] for i in range(n)]
while True:
NF = NF_generate(F, limit)
#节点分层
Vs = Vs_generate(NF, F, limit)
#构造辅助网络,改进分层网络（在传统定义上，只保留有后续边的节点，主要用于preflow_push）
flag, AN, Vs = AN_generate(Vs, NF)
# 若汇点不在分层网络中，则已得到最大流
if flag == -1:
break
if type_push == "DFS":
F = F_iteration_DFS(F, AN)
elif type_push == "preflow_push":
F = F_iteration_preflow_push(F, AN, Vs)
else:
F = F_iteration(F, AN)
return F
def Experiment(times):
time_Ford_Fulkerson_opt = []
time_A = []
time_IP = []
time_dinic = []
time_preflow_push = []
for nums in range(2,times):
limit = limit_generate(nums, 1, 100)
t0 = time.clock()
F_Ford_Fulkerson_opt = Ford_Fulkerson_opt(limit)
t1 = time.clock()
F_A, R, N = A_matrix(limit)
t2 = time.clock()
F_IP = IP(limit, 100)
t3 = time.clock()
F_dinic = dinic(limit, "DFS")
t4 = time.clock()
F_preflow_push = dinic(limit, "preflow_push")
t5 = time.clock()
time_Ford_Fulkerson_opt.append(t1 - t0)
time_A.append(t2 - t1)
time_IP.append(t3 - t2)
time_dinic.append(t4 - t3)
time_preflow_push.append(t5 - t4)
print("\ntime_Ford_Fulkerson_opt:", time_Ford_Fulkerson_opt)
print("\ntime_A:", time_A)
print("\ntime_IP:", time_IP)
print("\ntime_dinic:", time_dinic)
print("\ntime_preflow_push:", time_preflow_push)
if __name__ == "__main__":
# Experiment(20)
#点的个数
nums = 5
#流量限制，无向图
limit = limit_generate(nums, 1, 100)
LL_print("limit", limit)
# 改进的Ford-Fulkerson
t0 = time.clock()
F_Ford_Fulkerson_opt = Ford_Fulkerson_opt(limit)
#邻接矩阵方法
t1 = time.clock()
F_A, R, N = A_matrix(limit)
Route_print("Route", R, N)
print("\n最大流:" + str(sum(N)))
#整数规划方法：1邻接建模、2点边建模
t2 = time.clock()
F_IP = IP(limit, 100)
t3 = time.clock()
F_IP2 = IP2(limit, 100)
#dinic方法
t4 = time.clock()
F_dinic = dinic(limit, "DFS")
#preflow_push方法
t5 = time.clock()
F_preflow_push = dinic(limit, "preflow_push")
t6 = time.clock()
LL_print("limit", limit)
LL_print("F_Ford_Fulkerson_opt", F_Ford_Fulkerson_opt)
LL_print("F_A", F_A)
LL_print("F_IP", F_IP)
LL_print("F_IP2", F_IP2)
LL_print("F_dinic", F_dinic)
LL_print("F_preflow_push", F_preflow_push)
print("sum_F_Ford_Fulkerson_opt: ", sum([sum(F_Ford_Fulkerson_opt[i]) for i in range(nums)]))
print("sum_F_A: ", sum([sum(F_A[i]) for i in range(nums)]))
print("sum_F_IP: ", sum([sum(F_IP[i]) for i in range(nums)]))
print("sum_F_IP2: ", sum([sum(F_IP2[i]) for i in range(nums)]))
print("sum_F_dinic: ", sum([sum(F_dinic[i]) for i in range(nums)]))
print("sum_F_preflow_push: ", sum([sum(F_preflow_push[i]) for i in range(nums)]))
print("sum_F_Ford_Fulkerson_opt[0]: ", sum(F_Ford_Fulkerson_opt[0]))
print("sum_F_A[0]: ", sum(F_A[0]))
print("sum_F_IP[0]: ", sum(F_IP[0]))
print("sum_F_IP2[0]: ", sum(F_IP2[0]))
print("sum_F_dinic[0]: ", sum(F_dinic[0]))
print("sum_F_preflow_push[0]: ", sum(F_preflow_push[0]))
print("time_F_Ford_Fulkerson_opt: ", t1-t0)
print("time_F_A: ", t2-t1)
print("time_F_IP: ", t3-t2)
print("time_F_IP2: ", t4 - t3)
print("time_F_dinic: ", t5-t4)
print("time_F_preflow_push: ", t6 - t5)
扫一扫
6万+
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
