Python基础图论算法

本来下学期在学姐的强力安利之下选了算法这个课，但是取消了，于是在家打发时间看了edX上的一个法国人讲的算法网课。主要讲一些基础的图论算法，结合一个老鼠走迷宫的问题，用Python 写写程序。

image.png

个人感觉这个课对Python新手还是很友好的，因为他让你写的程序每一行都给你了注释。不想花钱搞verified，又觉得过段时间记录丢失有点亏，于是总结一下，留给自己看。
后来在回来的飞机上看了看5047 Intelligence System的PPT，发现这些算法全在在5047的Week 2 lecture里。我看了一个多星期的内容居然这个lecturer要用两个小时加50页ppt讲完。所以基本上只是提到了概念，没有做深入探讨 (毕竟这篇文章加上程序4000字)。

一， 先说算法

这里涉及到的算法有：DFS, BFS, Dijkstra, Brute force, Back tracking, Greedy。
正常的一个图，由边（Edge）和节点（Vertex）组成，可以用矩阵表示（Corresponding Matrix）。如果不考虑边的长度，两个点之间连接算1，不连接算0，这种叫adjacency matrix。这个矩阵横坐标代表从V1到Vn，纵坐标也是V1到Vn，关于主对角线对称（往返），并且对角线上全是0（自己到自己）。

image.png

Tree是在图的基础上得到的，要求cycle-free 。随便定义一个点为root，这样其他点就有了parent和children 的关系。因为cycle-free，一个children只能对应一个parent，这样从root到每一个节点的path就被唯一确定。有点像9135里网络的Routing Table。根据图我们可以通过BFS或者DFS得到图的Spinning Tree。这个过程叫做graph traversal ，这个结果不唯一。

image.png
DFS指的是deep first search，深度优先，就是一条路走到黑，没路了再回头。每连接一个节点，下一个节点选相邻的没连接的，如果相邻的节点都已经连接过了，就返回，到上一步里选一个相邻的没连接的，直到所有节点都连接到这个tree上。DFS产生的tree的特点是分支尽可能的少。
image.png
BFS指的是breadth first search，宽度优先，就是先把周围的连接上，再去连接更远的点。从一个节点开始，先连接和他相邻的节点，再连接隔一个hop的到和他相邻的，再连接隔两个hop的到隔一个hop的节点上，直到所有点都连接到这个tree上。BFS生成的tree有一个特点，root到任意一个节点经过的hop都是最少的情况。
image.png

如果每个边的长度都是1，那么BFS产生的tree就是从root到任意一个点的最短距离的情况。但是如果考虑比较一般的情况，边的长度是不一定相等的正数，这时候要找到一个点的最短距离就要用Dijkstra 算法。

image.png
先考虑最简单的情况，从v1到v3的距离。v2和 v3之间的距离如果大于2，显然最短距离是v1到v3，如果小于2，最短距离是v1-v2-v3。这个最短路径需要进一步的比较才能确定。但是，如果我们要找的是从v1到v2的最短距离，因为所有边的距离都是正数，很显然不管v2到v3的距离是多少，v1直接到v2一定是最短距离，而不是经过v3。所以我们可以得出一个结论，一个节点和他相邻的节点里距离最短的节点之间的路径，就是这两个节点之间的最短路径。
对于比较复杂的情况，使用迭代的思想。从root开始，找相邻节点中距离最短的，然后把已经确定的最短路径上的节点和root节点当作一个整体，更新这个大的root节点到相邻节点之间的距离为到实际root节点的距离，root整体节点不断扩大，不断迭代，直到包含所有的节点，就能找到每一个节点距离都最近的tree。
image.png
image.png
image.png
例子如上图所示，能够确定v1到v3一定是最短距离，于是把v1和v3当作一个整体，一个大的root节点，然后整体节点到相邻节点v2的最短距离从5更新为4。这时候包含v1和v3的大的root节点相邻的路径到v2是最短（到v2更新为4，到v7更新为5），下一个被包进来的节点就是v2，以此类推，实现迭代，最终找到root到每个节点的最短路径。
以上考虑的都是从一个节点到其他节点的距离问题。下面考虑一个更复杂的问题Traveling Salesman Problem，简称TSP。这个问题是在找从一个点开始，把图上的所有节点走一遍的最短距离。
image.png
image.png
这个问题图上如果有N个节点，那么就要列举计算(N-1)!种情况，这个数字非常大，所以这类问题没有一个有效的算法去解决。所谓Brute force就是计算出所有可能的路径，找出最小的。Back tracking在此基础上稍稍优化，记录一个当前找到的最短路径，在计算下一条路径的时候如果算到一半发现总的距离已经比这个最短路径大了，就不再进行这条路上后面的计算，返回，算下一条路。原理很简单，只是稍作优化。

Greedy算法就是说每一步都是局部最优（locally optimal），对整体最有结果进行逼近，不一定能找到整体的最优结果，但是高效。放到这个TSP问题上，就是我在v1的时候走到v4是最短距离，在v4的时候走到v3是最短距离，从而认为v1-v4-v3-v2是最短距离，不再进行其他判断。这个结果显然不是最优的，但是高效。其他问题上也可以使用greedy算法，比如要凑齐一定的钱用最少的硬币。

image.png

二，再来写程序

以上内容的原理其实很好理解，写程序就比较难了。
在Python里表示一个图，使用的是字典（Dict）的嵌套。先把图上点分好坐标，作为外层字典的key，每个点的value值包含内层字典，表示其他点到这个点的距离。

maze_graph = {
(0,0): {(0,1):1,(1,0):1},
(0,1): {(0,2):1,(0,0):1},
(1,0): {(1,1):1,(0,0):1},
(1,1): {(1,2):1,(1,0):1},
(0,2): {(0,1):1,(1,2):1},
(1,2): {(0,2):1,(1,1):1}
}

写DFS需要用到一个数据结构，LIFO Stack（last in first out）。最后放到这个stack里（Push）的元素最先出来（Pop）。这个东西虽然Python里有，可以直接用，但也可以用List实现。从root节点开始，每次先pop一个节点，当作走过的，再把相邻的没走过的节点push进这个stack，先push进去的会最后pop出来，这样同一级的节点就会最后走。直到所有的节点都被pop，整个图就按照DFS的方式走了一遍。

LIFO_list = list()
def LIFO_push(LIFO_list,element):
LIFO_list.append(element)
def LIFO_pop(LIFO_list):
return LIFO_list.pop(-1)

这里定义一个list，用来存已经走过的节点，并且定义方法来加加上已经走过的节点，判断一个点是否已经走过，即是否在这个list中。

def add_to_explored_vertices(explored_vertices,vertex):
       explored_vertices.append(vertex)

def is_explored(explored_vertices,vertex):
       return vertex in list(explored_vertices)

以下是实现这个DFS算法的方法，参数是图和root节点，返回值是已经走过的节点和每个节点的parent的字典。Python比Java在编程上的一个很明显的改进就是支持一个方法可以有很多个返回值，读取这些返回值时只要把变量按照顺序排好就可以。

def DFS(maze_graph, initial_vertex) :
# explored vertices list
      explored_vertices = list()
#LIFO stack
      queuing_structure = list()
#Parent Dictionary
      parent_dict = dict()
      LIFO_push(queuing_structure,(initial_vertex,None)) # push the initial vertex
     while len(queuing_structure) > 0: # while queuing_structure is not empty:
# current_vertex,parent = queuing_structure.pop()
# if the current vertex is not explored
# add current_vertex to explored vertices
# use parent_dict to map the parent of the current vertex
# for each neighbor of the current vertex in the maze graph:
# if neighbor is not explored:
# push the tuple (neighbor,current_vertex) to the queuing_structure
          current_vertex,parent =queuing_structure.pop()
          if not is_explored(explored_vertices,current_vertex):
          add_to_explored_vertices(explored_vertices,current_vertex)
          parent_dict.update({current_vertex:parent})
         for neighbor in maze_graph[current_vertex]:
                      if not is_explored(explored_vertices,neighbor):
                          LIFO_push(queuing_structure,(neighbor,current_vertex))
return explored_vertices,parent_dict

写BFS需要用到另外一个数据结构，FIFO Queue（first in first out），先放进去的元素后出来。这个和DFS的原理差不多，也是先pop出来一个节点，当作走过的，再把相邻的没走过的节点push进这个queue，因为先放进去的会先出来，这样就会先走同一级的节点，直到所有的节点都走遍。

def BFS(maze_graph, initial_vertex) :
# explored vertices list
        explored_vertices = list()
#LIFO stack
        queuing_structure = list()
#Parent Dictionary
       parent_dict = dict()
       FIFO_push(queuing_structure,(initial_vertex,None)) # push the initial vertex to the while      
       len(queuing_structure) > 0: # while queuing_structure is not empty:
# current_vertex,parent = queuing_structure.pop()
# if the current vertex is not explored
# add current_vertex to explored vertices
# use parent_dict to map the parent of the current vertex
# for each neighbor of the current vertex in the maze graph:
# if neighbor is not explored:
# push the tuple (neighbor,current_vertex) to the queuing_structure
       current_vertex,parent = FIFO_pop(queuing_structure)
       if current_vertex not in explored_vertices:
            explored_vertices.append(current_vertex)
            parent_dict.update({current_vertex:parent})
            FIFO_push(queuing_structure,(current_vertex,parent))
            for neighbor in maze_graph[current_vertex]:
     if neighbor not in explored_vertices:
                        FIFO_push(queuing_structure,(neighbor,current_vertex))
return explored_vertices,parent_dict

以下这个方法用来利用DFS或者BFS 的结果找到图上从一个点到任意一个点的路径。

def create_walk_from_parents(parent_dict,initial_vertex,target_vertex):
        relist = list()
        relist.append(target_vertex)
        while(initial_vertex != target_vertex):
    key = parent_dict[target_vertex]
    relist.append(key)
     target_vertex = key
     relist.reverse()
        if relist != []:
              del relist[0]
        return relist

initial_vertex = (0,0)
target_vertex = (0,2)
explored_vertices,parent_dict = DFS(maze_graph,initial_vertex)
route = create_walk_from_parents(parent_dict,initial_vertex,target_vertex)

在写Dijkstra 算法的时候，需要用到另一个数据结构，min-heap。Min-heap可以理解为一个字典，里面都是（key, value）一对一对的。这个min-heap主要用两个操作，add-or-replace和remove，add-or-replace就是在添加一个新的(key, value)的时候，如果这个key在min-heap里已经有了，value只有小于min-heap中的，才能更新value。这样可以保证value始终是最小值，并且key值不重复。因为Dijkstra需要找到最短路径，我们就可以把节点作为key，距离作为value，不断更新这个value，就能找到最短路径。
但是在做这个算法的时候，因为需要得到routing table，节点的parent也要加上。所以我们需要定义一个triplet，(vertex, weight, parent)，在比较value的时候，只比较weight。
先定义一个heap的pop方法，得到第一个元素。

# heap_pop function returns the first element of the list implementing the heap, providing the heap is not empty
def heap_pop(heap):
    if heap != []:
        vertex,weight,parent = heap.pop(0)
        return (vertex, weight, parent)
else:
       raise

因为要比较weight的值，即第二个元素，最后在sort中也要按照weight排序。定义一个方法，返回tuple中的第二个元素

def fun(s):
       return s[1]

之后定义这个add_or_replace方法。

def heap_add_or_replace(heap,triplet):
       i = 0
       while i < len(heap):
              if (heap[i])[0] == triplet[0]:
                    if (heap[i])[1] > triplet[1]:
                         heap[i] = triplet
                         heap.sort(key = fun)
                    return
              i = i+1
       heap.append(triplet)
       heap.sort(key = fun)

def is_explored(explored_vertices,vertex):
       return vertex in explored_vertices
def add_to_explored_vertices(explored_vertices,vertex):
       explored_vertices.append(vertex)
def Dijkstra(maze_graph,initial_vertex):
# Variable storing the exploredled vertices vertexes not to go there again
       explored_vertices = list()
# Stack of vertexes
       heap = list()
#Parent dictionary
       parent_dict = dict()
# Distances dictionary
       distances = dict()
# First call
       initial_vertex = (initial_vertex, 0, initial_vertex)#vertex to visit, distance from  
       heap_add_or_replace(heap,initial_vertex)
       while len(heap) > 0:
# get the triplet (vertex, distance, parent) with the smallest distance from heap # if the vertex of the triplet is not explored:
# map the vertex to its corresponding parent
# add vertex to explored vertices.
# set distance from inital_vertex to vertex
# for each unexplored neighbor i of the vertex, connected through an edge of # add (i, distance + wi, vertex) to the heap
#
                vertex,distance,parent = heap_pop(heap)
               if not is_explored(explored_vertices,vertex):
                     parent_dict.update({vertex:parent})
                     add_to_explored_vertices(explored_vertices,vertex)
                     distances[vertex] = distance
#distances.update({vertex:(distance+d)})
                     for neighbor in maze_graph[vertex]:
                             if neighbor not in explored_vertices:
                                   wi = (maze_graph[vertex])[neighbor]
                                   heap_add_or_replace(heap,(neighbor,wi + distance,vertex))

       return explored_vertices, parent_dict, distances

最后把图的尺寸作为变量，画一个计算时间随图的尺寸变量变化的关系图。

from imports import maze
import time
import math
%matplotlib inline
import matplotlib.pyplot as plt
maze_size=[]
time_to_execute=[]
maze_shape=[(1,2),(2,3),(3,3),(5,3),(6,7),(8,9),(11,10),(12,15),(16,18),(19,22),(21,23)]
for i in range(len(maze_shape)):
       length,width=maze_shape[i]
       maze_graph=maze.generate_maze(length,width,0,True,False,0.5,5,"",0)
       start=time.time()
       Dijkstra(maze_graph[3],(0,0))
       time_to_execute.append(time.time()-start)
       maze_size.append(length*width)
plt.plot(maze_size,time_to_execute,color="blue")

image.png

可以大致看出来程序运行时间基本和节点的数量是一个二次曲线的关系。这里涉及到算法复杂度的讨论。Dijkstra 算法的复杂度主要依赖与最小堆的实现方法。我们这个最小堆算每一个节点的距离的时候线性搜索了堆中的所有元素，所以复杂度是O(N^2). 经过对这个堆的算法优化，让他不是线性地去找最小值，可以进行优化，减少算法复杂度，实现O(NlogN)。这个比较难懂，不做探讨。
解决TSP问题的算法需要计算每一条路径，总路径的条数是(N-1)!，因此Brute force算法的复杂度就是O((N-1)!)。这个算法的复杂度大于N的多项式，不能在多项式时间内求解，这种问题就叫做NP问题（Non-deterministic Polynomial）。
Python里连接两个list可以用“+”，也可以用append()，区别是“+”不会改变原来的list()。

def create_vertices_meta_graph(piece_of_cheese, player_location):
#
       return piece_of_cheese + [player_location]
#

定义方法，把side的距离添加到adjacency矩阵里。

import utils
def create_edge_weight_maze_graph(maze_graph,vertices):
       adjacency_matrix={}
# for each initial_vertex in vertices:
# considere this vertex as source vertex
# use this source vertex and maze_graph to browse the graph with dijkstra algorithm
# for each vertex in vertices:
# use adjacency_matrix to store distances between source vertex and each vertex # remember to not store the distance from the source vertex to itself.
#
       for iv in vertices:
              source = iv
              explored_vertices,parent_dict,dist_dict = utils.Dijkstra(maze_graph,source)
              adjacency_matrix.update({source:{}})
              for v in vertices:
                    if v != source:
                         adjacency_matrix[source].update({v:dist_dict[v]})
       return adjacency_matrix

定义Brute force算法的方法，找到最短路径和最短距离。

def auxbf(current_walk,best_walk,adjacency_matrix,vertices,current_distance,best_distance):
# if length of current walk is larger than the order of the graph:
# if moreover the current distance is shorter than the best distance obtained so # update the value of the best distance
# update the corresponding best walk
# otherwise:
# for any vertex in the graph:
# if the vertex is not in the current walk:
# obtain potential values for best walk and best distance by recursively # if potential distance is shorter than best distance:
# update the value of the best distance
# update the corresponding best walk
#
       if len(current_walk) > len(adjacency_matrix)-1:
           if current_distance < best_distance:
                  best_distance = current_distance
                  best_walk = current_walk
        else:
               for v in adjacency_matrix:
                     if v not in current_walk:
                         walk_temp,dis_temp =   
                         auxbf(current_walk+[v],best_walk,adjacency_matrix,vertices,current_distance+   
                         adjacency_matrix[current_walk[-1]][v],best_distance)
                         if dis_temp < best_distance:
                             best_distance = dis_temp
                             best_walk = walk_temp

return best_walk,best_distance

调用上述三个方法进行brute force TSP问题的计算。

def bruteforceTSP(maze_graph,pieces_of_cheese,player_location):
# first we compute the vertices of the meta_graph:
       vertices=create_vertices_meta_graph(pieces_of_cheese, player_location)
# then we create the adjacency matrix of the meta graph
        adjacency_matrix = create_edge_weight_maze_graph(maze_graph,vertices)
# now we can start defining our variables
# current_distance is the length of the walk for the current exploration branch
        current_distance = 0
# current_walk is a container for the current exploration branch
         current_walk = [player_location]
# best_distance is an indicator of the shortest walk found so far
         best_distance = float('inf')
# best_walk is a container for the corresponding walk
        best_walk = []
# we start the exploration:
       best_walk, best_distance =  auxbf(current_walk,best_walk,adjacency_matrix,
       pieces_of_cheese,current_distance,best_distance)
       return best_walk, best_distance

在brute force的基础上稍加改动，加入一个和当前最短距离的判断，就能定义back tracking 方法

def auxbt(current_walk,best_walk,adjacency_matrix,vertices,current_distance,best_distance):
# First we test if the current walk have gone through all vertices
# if that is the case, we compare the current distance to the best distance
# and in the case it is better we update the best distance and the best walk
# if the current_walk is not finished we see if the current distance is lower than # if that is the case, for each possible vertex not explored,
# we add it and call ourself recursively
##
YOUR CODE HERE
#
if (len(current_walk) == len(vertices)):
    if (current_distance < best_distance):
        best_distance=current_distance
      best_walk=current_walk
else:
        if current_distance < best_distance:
            for v in adjacency_matrix:
                  if v not in current_walk:
                      walk_temp,dis_temp =  auxbt(current_walk+[v],best_walk,adjacency_matrix,
                      adjacency_matrix[current_walk[-1]][v],best_distance)
                      if dis_temp < best_distance:
                          best_distance = dis_temp
                          best_walk = walk_temp
return best_walk,best_distance

用time 计算两个算法对不同尺寸的图的计算时间，并绘图比较。

import time
%matplotlib inline
import matplotlib.pyplot as plt
maze_shape = [(3,3),(5,6),(6,7),(9,10),(11,13)]
poc = [2,3,4,5,6]
bruteforce_time = list()
backtracking_time = list()
maze_size=list()
for i in range(len(maze_shape)):
      width,height = maze_shape[i]
      _,_,_,maze_graph = maze.generate_maze(width,height,0,True,False,0.5,5,"",0)
      pieces_of_cheese,player1_location,_ = maze.generate_pieces_of_cheese(
      poc[i], width, height, False, None, None, False)
      start = time.time()
      bruteforceTSP(maze_graph,pieces_of_cheese,player1_location)
      bruteforce_time.append(time.time()-start)
      start = time.time()
      backtrackingTSP(maze_graph,pieces_of_cheese,player1_location)
      backtracking_time.append(time.time() - start)
      maze_size.append(width*height)
plt.plot(maze_size,bruteforce_time,color="blue")
plt.plot(maze_size,backtracking_time,color="green")

分别是3x3的图里2个目标点，5x6的图里3个，6x7的图里4个，9x10的图里5个，11x13的图里6个的情况。可以看出back tracking在brute force的基础上有所优化，尤其是目标点多了之后。

image.png

最后是greedy算法，这个背景是老鼠要找迷宫中的所有奶酪，就每次都去找最近的，实现对总距离最短的近似实现。每一次的最近距离使用Dijkstra 算法找。

def choose_target(distances,pieces_of_cheese):
    distance = float('inf')
    next_target = None
    # for each vertex in pieces of cheese
    #    compare its distance from player location with distance variable value.
    #    update distance variable with the smallest value
    #    keep the nearest vertex from player location
    # returns the nearest vertex to be next target
    for vertex in pieces_of_cheese:
        if distances[vertex] < distance or distance is float('inf') :
            distance = distances[vertex]
            next_target = vertex    
return next_target

import utils
movements = list()
def movements_greedy_algorithm(maze_graph,pieces_of_cheese,player_location):
        # get distances using Dijkstra algorithm
    # get next_target using choose_target function
    # use A_to_B function to get a list of movements that should be done to reach the nearest piece of cheese
    # from player location
    explored_vertices,parent_dict,distance_dict = utils.Dijkstra(maze_graph,player_location)
    next_target = choose_target(distance_dict,pieces_of_cheese)
    movements = utils.A_to_B(maze_graph,player_location,next_target,parent_dict)
return movements,next_target

最后，画一个程序执行时间和图的节点数的关系图。

image.png