力引导图python实现 force directed layout

force directed layout 力引导图python实现

---

实验内容：

• force directed layout编程实现
• 探讨三个force directed layout算法的加速策略

1. A brief introduction to force directed layout

1.1 维基百科

Force-directed graph drawing algorithms are a class of algorithms for drawing graphs in an aesthetically-pleasing way. Their purpose is to position the nodes of a graph in two-dimensional or three-dimensional space so that all the edges are of more or less equal length and there are as few crossing edges as possible, by assigning forces among the set of edges and the set of nodes, based on their relative positions, and then using these forces either to simulate the motion of the edges and nodes or to minimize their energy.

While graph drawing can be a difficult problem, force-directed algorithms, being physical simulations, usually require no special knowledge about graph theory such as planarity.

1.2 理论基础

力导向图利用模拟物理世界粒子间相互作用力的思想，以下边两个理论为基础：

• 粒子间库仑力：
  $$F=k_r\frac{q_1{q}_2}{l^2}$$
  在实际应用中，我们认为$q_1=q_2=1$，$l$是两个点之间的距离，$k_r$是人为设定的超参数。在算法中所有点之间都要计算排斥力。

• 弹簧胡克定律：
  $$F=k_s(l-r)$$
  模拟弹簧力，其中$r$是弹簧自然长度，$k_s$是弹簧的劲度系数，$l$是两点之间距离。在算法中只有邻接的两个点之间才计算弹簧力。

1.3 优势

• **布局结果质量比较高：**对于点数在50-500之间的图布局结果较好
• **灵活性：**力导向算法可以轻松进行调整和扩展，以满足其他审美标准。这使它们成为最通用的图形绘制算法类
• **算法行为容易预测：**因为理论基础来源于物理
• **简单：**经典的力导向图代码量不大
• **互动性：**用户可以了解算法的执行过程，绘制每一步迭代的结果图

1.4 缺点

• 时间复杂度高
• 可能形成局部最优

2. python实现效果

def init():    
    for i in range(0,Node_num):#随机生成点坐标,初始化力
        posx=random.uniform(0,original_max_posx)
        posy=random.uniform(0,original_max_posy)#初始化点坐标和受力
        Node_position[i]=(posx,posy)
        Node_force[i]=(0,0)
    for i in range(0,int(Node_num/2)-1):#随机生成边
        index=random.randint(0,i)
        Edge.append((i,index))
        index=random.randint(0,i)
        Edge.append((i,index))
    for i in range(int(Node_num/2),Node_num):#随机生成边
        index=random.randint(int(Node_num/2),Node_num-1)
        Edge.append((i,index))

• 生成随机生成无向无环连通图：

  在下图中，节点的大小用其degree表示，度越大，节点的size越大。随机生成初始位置，经过迭代后结果如下。可以看到最后degree较大的节点呈现中心收敛现象。参数设置为：$k_r=6,k_s=0.3$，经过多次尝试，这组参数效果最好。一般而言，库仑力系数$k_r$要比$K_s$大得多，因为在后续迭代过程中可以大概保持图的全局性，不会收缩成一团。

  

• 生成随机生成无向无环非连通图：

  ​ 在实验中发现一个有趣的现象，force direction layout可以展现图的层次结构，下边是几组实验结果

  球状：
  

  树状：
  

  惊奇的发现force direction layout可以把不同连通分支分开，并且叶子节点大多分布在边缘，这可以用库仑力解释：

  

  也可以看到明显的层次结构：

  

  我们加上颜色通道来表现连通分支：

  

  对于有环图，经过一定次数的迭代最后收敛到很好的结果：

  

3. 加速策略讨论

3.1 模拟退火法

​ 模拟退火算法来源于固体退火原理，即高温状态的固体渐渐冷却，固体内部粒子动能逐渐变小，最终在常温状态下达到平衡，此时动能最小。在模拟退火算法有一个温度的概念，“温度”从一个初值逐渐减小为0，与此同时限制节点的最大偏移量也逐渐减小，最终达到在恰当时机控制算法终止的目的。

​ 在位置更新这一函数体内，让每一次迭代粒子运动的最大值减小。计算整个图粒子位移总和，这样。允许移动的最大步长设置为关于迭代次数的减函数，这里简单的设置为$y=\frac{1}{x}$。那么在初期节点允许移动的最大距离很大，后期很小，有很大概率停下来。

• 对节点的偏移总量进行连续采样，当图布局变化不大时结束迭代。

  

• 结果：

  看200个点的迭代：
  

  

• 200个点无退火优化，看起来虽然比有优化的分层次更明显，但是上述优化的结果也可以接受，最重要的是退火法减少了用时。
  

3.2 合并同一连通分支内的紧凑节点

​ 对于多连通分支图布局来说可以通过合并同一连通子图中的节点为一个节点的思想减少计算节点之间的斥力。**比如下图的红色连通分支认为已经收敛，那么我们可以用它的坐标中心作为其代表点的坐标。**这样点的个数就变少了，计算的次数也相应减少。

​ 进一步的，我们为节点的位置变化设置一个阈值，如果连续n次迭代节点位置变化的方差小于这个阈值（可以理解为节点的位置不变化了，稳定了）。那么我们可以将此类节点合并成新的节点。这一想法应该对于区分不同连通分量比较有效。

3.3 Barnes–Hut simulation

​ 对于节点数量较大的图布局时退火法运行时间仍然很大。我们可以考虑新的方法：

​ Barneshut算法主要思想是根据节点位置距离建立一颗树，将邻近多个的节点看成一个超级节点，从而减少了斥力计算环节中需要计算任意两点间的斥力的复杂度，这部分的时间复杂度由$O(n^2)$降为了$O(nlogn)$，最终引入Barneshut算法后使得力导向算法的时间复杂度降为$O(k(nlogn+m))$。

​ Barneshut算法是很巧妙的方法，将邻近区域的节点分组合并，广泛用于n-body仿真。它递归地将节点集合存储在四叉树结构中。顶点代表整个区域，它的质量为所有节点质量之和，它的位置为所有子节点的质心位置。这个算法之所以快是因为我们不需要去计算每一个组body里面的节点。

​ Barneshut算法步骤：

• 创建根节点body，不断地按节点位置将所有节点划分到body的四个象限，由此建立树结构。

• 遍历计算每一个节点与树结构之间的斥力，若当前节点的位置与树结构body节点的质心位置足够远（s/d<0.5),则将所有的作用力施加于根节点上；若不足够远，则递归地计算当前节点与body节点的子节点的斥力。
  

4. python 代码

"""
力引导图布局
"""
import random
import math
from matplotlib import pyplot as plt
import networkx as nx
import time
import numpy as np
#模型参数
K_r = 6
K_s = 0.3
L = 5
delta_t = 50
MaxLength = 30
iterations=200
color=['red','green','blue','orange'] 
Node_num=200
original_max_posx=40
original_max_posy=40
#图形容器
Edge=[]
Node_force={}
Node_position={}
#节点大小用来反映节点的度
Node_degree=[]

#与采样退火有关的参数
Displacement_list=[]#用于采样的列表
scale=3#采样的范围；

def init():    
    for i in range(0,Node_num):#随机生成点坐标,初始化力
        posx=random.uniform(0,original_max_posx)
        posy=random.uniform(0,original_max_posy)#初始化点坐标和受力
        Node_position[i]=(posx,posy)
        Node_force[i]=(0,0)
    for i in range(0,int(Node_num/2)-1):#随机生成边
        index=random.randint(0,i)
        Edge.append((i,index))
        index=random.randint(0,i)
        Edge.append((i,index))
    for i in range(int(Node_num/2),Node_num):#随机生成边
        index=random.randint(int(Node_num/2),Node_num-1)
        Edge.append((i,index))

def compute_repulsion():#计算每两个点之间的斥力
    for i in range(0,Node_num):
        for j in range(i+1,Node_num):
            dx=Node_position[j][0]-Node_position[i][0]
            dy=Node_position[j][1]-Node_position[i][1]
            if dx!=0 or dy!=0:
                distanceSquared=dx*dx+dy*dy
                distance=math.sqrt(distanceSquared)
                R_force=K_r/distanceSquared
                fx=R_force*dx/distance
                fy=R_force*dy/distance#更新受力
                fi_x=Node_force[i][0]
                fi_y=Node_force[i][1]
                Node_force[i]=(fi_x-fx,fi_y-fy)
                fj_x=Node_force[j][0]
                fj_y=Node_force[j][1]
                Node_force[j]=(fj_x+fx,fj_y+fy)

def compute_string():
    for i in range(0,Node_num):#取出其邻居
        neighbors=[n for n in G[i]]#对每一个邻居，计算斥力；j
        for j in neighbors:
            if i < j:
                dx=Node_position[j][0]-Node_position[i][0]
                dy=Node_position[j][1]-Node_position[i][1]
                if dx!=0 or dy!=0:
                    distance=math.sqrt(dx*dx+dy*dy)
                    S_force=K_s*(distance-L)
                    fx=S_force*dx/distance
                    fy=S_force*dy/distance#更新受力
                    fi_x=Node_force[i][0]
                    fi_y=Node_force[i][1]
                    Node_force[i]=(fi_x+fx,fi_y+fy)
                    fj_x=Node_force[j][0]
                    fj_y=Node_force[j][1]
                    Node_force[j]=(fj_x-fx,fj_y-fy)
                
def update_position(times):#更新坐标
    Displacement_sum=0
    for i in range(0,Node_num):
        dx = delta_t*Node_force[i][0]
        dy = delta_t*Node_force[i][1]
        displacementSquard = dx*dx + dy*dy
        #随迭代次数增加，MaxLength逐渐减小；

        current_MaxLength = MaxLength/(times+0.1)

        if( displacementSquard >current_MaxLength):
            s=math.sqrt(current_MaxLength/displacementSquard)
            dx=dx*s
            dy=dy*s
        (newx,newy) = (Node_position[i][0]+dx, Node_position[i][1]+dy)
        Displacement_sum += math.sqrt(dx*dx + dy*dy) 
        Node_position[i]=(newx,newy)
    return Displacement_sum

if __name__ == '__main__':
    G=nx.Graph()
    G.add_nodes_from(list(range(0,Node_num)))
    init()#初始化节点得坐标和图中的边
    G.add_edges_from(Edge)
    #获得节点的度
    for i in range(0,Node_num):
        Node_degree.append(pow(G.degree(i),2))
    #获得联通子图
    connected_num=0
    connected_subgraph=[]
    for c in nx.connected_components(G):
        connected_num += 1
        nodeSet = G.subgraph(c)
        connected_subgraph.append(nodeSet)
    
    #原图
    nx.draw_networkx_nodes(G, pos = Node_position, node_size=20, node_color = 'red',alpha=0.8)
    nx.draw_networkx_edges(G, pos = Node_position,edge_color='lightblue',width=0.5) 
    plt.show()

    fig=plt.figure('不同迭代次数force-direction_layout')
    start =time.perf_counter()
    iteration_time=0
    for times in range(0,1+iterations):
        for i in range(0,Node_num):
            Node_force[i]=(0,0)
        compute_repulsion()
        compute_string()
        #记录本次迭代移动距离：
        Displacement_sum = update_position(times)
        Displacement_list.append(Displacement_sum)
        print(Displacement_sum)
        if len(Displacement_list)>scale:
            last = np.mean(Displacement_list[times-4:times-1])
            now = np.mean(Displacement_list[times-3:times])
            if (last-now)/last < 0.01:
                break
        iteration_time=times
    end = time.perf_counter()

    print('Running time: %s Seconds'%(end-start))
    print('最终迭代次数:',iteration_time)
    index=0
    for subgrap in connected_subgraph:
        sub_position = dict({i for i in Node_position.items() if i[0] in list(subgrap)})
        sub_degree = [degree+5 for (point,degree) in enumerate(Node_degree) if point in list(subgrap) ] 
        nx.draw_networkx_nodes(subgrap, pos = sub_position, node_size=sub_degree, node_color = color[index%4],alpha=0.8)
        nx.draw_networkx_edges(subgrap, pos = sub_position,edge_color='lightblue',width=0.5) 
        index+=1
    #受力得初始化：
    plt.show()