我做KDD比赛时学习KM算法的心路历程

import scipy.sparse as sp
import cvxpy as cp
import numpy as np

# 假设我们现在有司机和订单的id列表，并且有一些司机和订单匹配估算得到的价格（注意：并不是每个订单和司机都会匹配，有可能有些是没有匹配的，所以这个矩阵是一个稀疏的矩阵，因此要引入sicpy里的sparse来生成稀疏矩阵

row_num = 100  # 司机个数

col_num = 200  # 订单个数

val_num = 1000 # 权重矩阵里有值的个数，也就是初始提供的司机订单配对数

values = 1 + np.random.randn(val_num) * 2  # 随机生成值

rows = np.random.randint(0, row_num, val_num)  # 随机生成初始配对值所对应的横坐标

cols = np.random.randint(0, col_num, val_num)  # 随机生成初始配对值所对应的纵坐标

A_sp = sp.coo_matrix((values, (rows, cols)), shape=(row_num, col_num))  # 生成稀疏矩阵

# cvxpy求解凸优化问题
X = cp.Variable((row_num, col_num), boolean=True)  # 生成一个只有0和1的变量，维数和A_sp相同

obj = cp.Maximize(cp.sum(cp.multiply(A_sp, X)))  # 要优化的目标

constraints = [X @ np.ones((col_num, 1)) <=1, 
               np.ones((1, row_num)) @ X <=1]   # 满足的约束，右边的1其实起到的是向量的作用

problem = cp.Problem(obj, constraints)     # 生成problem

problem.solve()  # 求解

problem.value    # 输出最大值

probelm.X        # 输出对应的X的值

import scipy.sparse as sp
import cvxpy as cp
import numpy as np

# 假设我们现在有司机和订单的id列表，并且有一些司机和订单匹配估算得到的价格（注意：并不是每个订单和司机都会匹配，有可能有些是没有匹配的，所以这个矩阵是一个稀疏的矩阵，因此要引入sicpy里的sparse来生成稀疏矩阵

row_num = 100  # 司机个数

col_num = 200  # 订单个数

val_num = 1000 # 权重矩阵里有值的个数，也就是初始提供的司机订单配对数

values = 1 + np.random.randn(val_num) * 2  # 随机生成值

rows = np.random.randint(0, row_num, val_num)  # 随机生成初始配对值所对应的横坐标
cols = np.random.randint(0, col_num, val_num)  # 随机生成初始配对值所对应的纵坐标


# cvxpy求解凸优化问题

order_mat_row, order_mat_col = [None] * val_num, [None] * val_num
driver_mat_row, driver_mat_col = [None] * val_num, [None] * val_num

for idx, (did, oid) in enumerate(zip(rows, cols)):
    order_mat_row[idx], order_mat_col[idx] = oid, idx
    driver_mat_row[idx], driver_mat_col[idx] = did, idx

driver_mat = sp.coo_matrix((np.ones(val_num), (driver_mat_row, driver_mat_col)), shape=(row_num, val_num))

order_mat = sp.coo_matrix((np.ones(val_num), (order_mat_row, order_mat_col)), shape=(col_num, val_num))

# 上面这一段相当于建立分别建立司机和订单的稀疏矩阵，以司机为例，这个矩阵行是司机的id，列是所有司机订单的配对值，对于每一行来说，记录的是该司机在哪些配对值中出现过，出现为1，否则为空；订单矩阵的含义类似

# cvxpy求解凸优化问题

X = cp.Variable(val_num, boolean=True)  # 由原来的选哪些司机和订单变成直接是选哪些值

obj = values@X
constr = [
    order_mat @ X <= 1,    # 意味着每个订单只能匹配不超过一次
    
    driver_mat @ X <= 1,   # 每个司机只能匹配不超过一次
]
prob = cp.Problem(cp.Maximize(obj), constr)
prob.solve(solver=cp.GLPK_MI, glpk={'msg_lev': 'GLP_MSG_OFF', 'presolve': 'GLP_ON'})  # 可以暂时忽略这些参数

prob.value
prob.X

graph = {'A': set(['B', 'C']),
         'B': set(['A', 'D', 'E']),
         'C': set(['A', 'F']),
         'D': set(['B']),
         'E': set(['B', 'F']),
         'F': set(['C', 'E'])}

def dfs(graph, start):
    visited, stack = set(), [start]
    while stack:
        vertex = stack.pop()  # 默认从最后删除
        
        if vertex not in visited:
            visited.add(vertex)
            stack.extend(graph[vertex] - visited)   # 这里一定注意要除掉已经访问的顶点，不然会进入死循环
    return visited

dfs(graph, 'A') # {'E', 'D', 'F', 'A', 'C', 'B'}

def dfs(graph, start, visited=None):
    if visited is None:
        visited = set()
    visited.add(start)
    for next in graph[start] - visited:
        dfs(graph, next, visited)
    return visited

dfs(graph, 'C') # {'E', 'D', 'F', 'A', 'C', 'B'}

def dfs_paths(graph, start, goal):
    stack = [(start, [start])]
    while stack:
        (vertex, path) = stack.pop()
        for next in graph[vertex] - set(path):
            if next == goal:
                yield path + [next]   # 返回一个生成器
                
            else:
                stack.append((next, path+[next]))
                
list(dfs_paths(graph, 'A', 'F'))    # [['A', 'C', 'F'], ['A', 'B', 'E', 'F']]

def dfs_paths(graph, start, goal, path=None):
    if path is None:
        path = [start]
    if start == goal:
        yield path
    for next in graph[start] - set(path):
        yield from dfs_paths(graph, next, goal, path+[next])
        
list(dfs_paths(graph, 'C', 'F')) # [['C', 'F'], ['C', 'A', 'B', 'E', 'F']]

def bfs(graph, start):
    visited, queue = set(), [start]
    while queue:
        vertex = queue.pop(0)  # 表示删除第一个
        
        if vertex not in visited:
               visited.add(vertex)
               queue.extend(graph[vertex]-visited)
    return visited

bfs(graph, 'A') # {'B', 'C', 'A', 'F', 'D', 'E'}

from collections import deque

def bfs(graph, start):
    visited= set()
    queue = deque([start])
    while queue:
        vertex = queue.popleft()  # 表示删除第一个
        
        if vertex not in visited:
               visited.add(vertex)
               queue.extend(graph[vertex]-visited)
    return visited

bfs(graph, 'A') # {'B', 'C', 'A', 'F', 'D', 'E'}

def bfs_paths(graph, start, goal):
    queue = [(start, [start])]
    while queue:
        (vertex, path) = queue.pop(0)
        for next in graph[vertex] - set(path):
            if next == goal:
                yield path + [next]
            else:
                queue.append((next, path + [next]))
                             
list(bfs_paths(graph, 'A', 'F'))   # [['A', 'C', 'F'], ['A', 'B', 'E', 'F']] 

def shortest_path(graph, start, goal):
    try:
        return next(bfs_paths(graph, start, goal))
    except StopIteration:
        return None
    
shortest_path(graph, 'A', 'F') # ['A', 'C', 'F']

# 假设我们有一个n*m的profit矩阵，n<=m，实际操作的时候，永远把较小的作为行就可以了，当然你也可以在函数里通过判断行列大小来写一个变换的操作，这样就不用考虑矩阵的形状了，但是这里就不再增加工作量了，默认在函数外我就已经变换好了

import numpy as np

def km_dfs(profit): 
    
    nrows, ncols = profit.shape
    # 初始化顶标
    
    left_top_idx = np.max(profit, axis=1).tolist()   # 左边就是上面例子中的上方，取每行最大值
    
    right_top_idx = [0] * ncols                      # 右边就是上面例子中的下方，赋值0
    
    left_match = [None] * nrows        # 用来记录每个行标对应的匹配列标
    
    right_match = [None] * ncols       # 用来记录每个列标对应的匹配行标，如果列大于行，会留有None
    
    
    def find_aug_path_dfs(i):
        '''这里寻找增广路径是对于每个行标而言的，因为我们是要遍历行标，所以每个行标都要走一遍这个搜索流程，
        
        并且用了递归来维护堆栈数据，回想之前提到过的搜索连接节点的问题，我们只添加未搜索过的节点，不然会进入死循环
        '''
        
        left_visit[i] = True
        for j, weight in enumerate(profit[i]):
            if right_visit[j]: continue # 已被访问（解决递归中的冲突）
                
            # RecursionError: maximum recursion depth exceeded while calling a Python object
            
            gap = left_top_idx[i] + right_top_idx[j] - profit[i][j]  # 计算gap是用来判断哪个Y是在相等子图里的，并且记录非0的gap值用于更新顶标
            
            if gap == 0:  # abs(gap) < zero_threshold
                
                right_visit[j] = True
                
                if right_match[j] == None or find_aug_path_dfs(right_match[j]): 
                # j未被匹配，或虽然j已被匹配，但是j的已匹配对象有其他可选项
                
                    left_match[i] = j                    
                    right_match[j] = i

                    return True
            else:
                
                if gap < slack[j]: slack[j] = gap   # 通过这一步将递归中的slack值不断更新，保持对于每个列标来说我都是取了最小的那个（虽然可能不是同一个i对应的j）

        return False
    

    for i in range(nrows):   # 遍历每一个行标
        
        while True:
            slack = [np.inf] * ncols    # 初始化slack用于记录gap值
            
            # 初始化搜索指标，每搜索一个行标都要重新初始化visited集合，在更新完顶标之后，即使是同一个行标，也是要重新初始化visited集合的，这个逻辑要理清
            
            left_visit = [False] * nrows
            right_visit = [False] * ncols
            
            if find_aug_path_dfs(i): break  # 如果找到增广路径，则跳出while，进行下一个行标的搜索，但是如果没找到增广路径，即返回的是False，则进行下面的更新顶标操作，并且依然在while循环里，为该行标寻找增广路径，如果依然是False，再修改，直到返回True
                
            d = np.inf
            for j, slk in enumerate(slack):
              if not right_visit[j] and slk < d:
                d = slk         # 考虑没有被搜索过的Y当中最小的gap，因为只要是调用find_aug_path_dfs()函数，那么该行标X就会被标记为访问过
                
            for r in range(nrows):
                if left_visit[r]: left_top_idx[r] -= d
            for c in range(ncols):
                if right_visit[c]: right_top_idx[c] += d
       
    value = 0
    for i, j in enumerate(left_match):
        value += profit[i][j]
        
    return left_match, value


profit = np.array([[3, 4, 6, 4, 9], 
                   [6, 4, 5, 3, 8],
                   [7, 5, 3, 4, 2],
                   [6, 3, 2, 2, 5],
                   [8, 4, 5, 4, 7]])

km_dfs(profit)
# ([4, 2, 3, 1, 0], 29)

def km_bfs(profit):
    nrows, ncols = profit.shape
    # 初始化顶标
    
    left_top_idx = np.max(profit, axis=1).tolist()
    right_top_idx = [0] * ncols
    
    left_match = [None] * nrows
    right_match = [None] * ncols

    def find_aug_path_bfs(i):
        left_visit[i] = True
        queue = []  # 记录那些对于第i个节点来说所有可能的匹配节点（即所有有边的节点）
        
        for j, weight in enumerate(profit[i]):
            gap = left_top_idx[i] + right_top_idx[j] - profit[i][j]
            if right_visit[j]: continue  # 防止无穷递归
                
            if gap == 0:
                right_visit[j] = True
                queue.append(j)
                if right_match[j] == None:   # 如果j节点没有被匹配过，那么配对成功，否则就要往下走，继续搜索对于i来说其他可能的节点
                    
                    right_match[j] = i
                    left_match[i] = j
#                    print(i, j, 'success')

                    return True          
            else:
                slack.append(gap)
        for j in queue:                  # 这里就是搜索对i来说可能的节点，如果能搜到，则完成匹配，不能的话就说明找不到增广路径，需要修改顶标
            
            if find_aug_path_bfs(right_match[j]):
                right_match[j] = i
                left_match[i] = j
#                print(i, j, 'success')

                return True
        
#        print(i, 'fail')

        return False

    for i in range(nrows):
        while True:
            # 初始化搜索指标
            
            left_visit = [False] * nrows
            right_visit = [False] * ncols 
            slack = []
        
            if find_aug_path_bfs(i): break
            d = min(slack)
            for r in range(nrows):
                if left_visit[r]: left_top_idx[r] -= d
            for c in range(ncols):
                if right_visit[c]: right_top_idx[c] += d

    value = 0
    for i, j in enumerate(left_match):
        value += profit[i][j]
    
    return left_match, value


profit = np.array([[3, 4, 6, 4, 9], 
                   [6, 4, 5, 3, 8],
                   [7, 5, 3, 4, 2],
                   [6, 3, 2, 2, 5],
                   [8, 4, 5, 4, 7]])

km_bfs(profit)
# ([4, 2, 1, 0, 3], 29)

# from collections import deque   # 使用deque的话就把相应的注释给去掉，同时把上一行给注释掉

def km_bfs(profit):
    
    nrows, ncols = profit.shape
    # 初始化顶标
    
    left_top_idx = np.max(profit, axis=1).tolist()
    right_top_idx = [0] * ncols
    
    left_match = [None] * nrows
    right_match = [None] * ncols
    
    
    
    def find_aug_path_bfs(s):
        
        slack = [np.inf] * ncols
        pre = [None] * ncols   # 用来记录slack是对应于哪个一个x的匹配
        
        # 初始化搜索指标
        
        left_visit = [None] * nrows
        right_visit = [None] * ncols
        
        queue = [s] # 建立队列，用来记录当次没有被匹配上的y上一次对应的x，将要搜寻的x作为初始值
        
#        queue = deque([s]) 

        left_visit[s] = True
        
        
        def check(j):
            right_visit[j] = True
            
            if right_match[j] != None:
                queue.append(right_match[j])
                left_visit[right_match[j]] = True
                return False
            
            while j != None:
                right_match[j] = pre[j]
                tmp = left_match[pre[j]]
                left_match[pre[j]] = j
                j = tmp
            return True
            
            
        while True:
            while queue:
                i = queue.pop(0)     # i = queue.popleft()
                
                for j in range(ncols):
                    gap = left_top_idx[i] + right_top_idx[j] - profit[i][j]
                    if not right_visit[j] and slack[j] >= gap:
                        pre[j] = i
                        if gap != 0:
                            slack[j] = gap
                        elif check(j):
                            return
        
            delta = np.inf
            for j in range(ncols):
                if not right_visit[j] and delta >= slack[j]:
                    delta = slack[j]
            
            for j in range(ncols):
                if right_visit[j]:
                    right_top_idx[j] += delta
                else:
                    slack[j] -= delta
            
            for i in range(nrows):
                if left_visit[i]:
                    left_top_idx[i] -= delta
            
            for j in range(ncols):
                if not right_visit[j] and slack[j] == 0 and check(j):
                    return
          
    for i in range(nrows):
        find_aug_path_bfs(i)
    
    value = 0
    for i, j in enumerate(left_match):
        value += profit[i][j]
        
    return left_match, value


profit = np.array([[3, 4, 6, 4, 9], 
                   [6, 4, 5, 3, 8],
                   [7, 5, 3, 4, 2],
                   [6, 3, 2, 2, 5],
                   [8, 4, 5, 4, 7]])

km_bfs(profit)
# ([2, 4, 3, 1, 0], 29)