在数据可视化的过程中，绘制网络拓扑图是很重要的，它能清晰呈现一个复杂网络的结构，节点的重要性和关系。比如下面几张图：
下面这张图是我的软件绘制的：
这些都有一个共同的问题，就是如何让图绘制的更加美观？ 复杂的图结构，已经没法人工布局了。所以计算机自动布局，就成了一个有趣而重要的话题。我们将其称为布点算法。
基本原理
一个好的布点算法应该能尽量满足以下四个特点：
- 对称性：绘制网络对应将具有相同结构的子图围绕绘图中心平衡布局
- 连接角最大化原则，使同一节点任意两条边形成的角度尽量大
边交叉数量最小原则（最重要）：绘图时应尽量减少相互交叉边的数量。 直线边原则：边尽量直线，避免曲线，同时线尽可能短。其中，尽量避免交叉是最重要的，然而，并不是所有的图结构都能满足这个要求。同时，在性能上，还应当满足如下需求：
算法应当能处理尽可能复杂和庞大的图，比如上万甚至上百万的节点
- 运算速度还要尽量快，使其能够在短时间内完成操作
- 可交互性，当用户调整了某一个点的位置，或改变了一些参数，算法是否能够动态调节？
如何表示一个网络
图分为有向图和无向图，有权图和无权图。当图有权时，权重反映了节点间的关系，通常值越大，节点关系越紧密。
- 二维矩阵 最简单也最方便的表示方法，当网络结构巨大时，存储空间浪费是惊人的，而且如果图是无向图，则矩阵是对称的。但存取效率最高
- 临接表， 即将边表达为一个数组，
(x,y,dist) *N对于稀疏图来说，存储效率很高，但读写效率低。一种简单的优化方法是对x和y 分别进行排序，排序后按二分查找，即可实现快速存取。 - 实时计算 当网络图巨大，而计算关系本身的运算并不复杂时，可以采用，此时dist[x][y] 被转换为一个函数 GetDistance(x,y), 实时完成计算。这是一种很有趣且比较实用的方法，尤其是当网络图稀疏时。
为了简化，我们采用二维矩阵作为存取结构，相信读者可以很容易的将其改写为其他形式。
算法介绍
力导引算法
力导引布局的方法可以产生相当优美的网络布局，并充分展现网络的整体结构及其自同构特征。该方法最早由Eades在1984年提出。
基本思想是将网络看成一个顶点为钢环，边为弹簧的物理系统，系统被赋予某个初始状态以后，弹簧弹力（引力和斥力）的作用会导致钢环移动，这种运动直到系统总能量减少到最小值停止。每次循环都要计算任意一对点的斥力和相邻两个点的引力，计算复杂度为O（N2）。
基本上绝大多数算法都遵循着这样的原则，即：
- 将网络看成一个顶点为钢环，边为弹簧的物理系统
- 不断迭代，使整个系统的总能量达到最小
为了简化，弹簧的受力并不遵循胡克定律，而是自己定义的受力公式，同时决定引力只在相邻节点间。因此，最基本的弹簧算法的伪代码如下：
随机分布G的节点
重复m次
{
计算作用在各个节点上的力
按照力移动节点的位置
}
根据节点位置绘图
衡量一个布点算法的美学标准，可以用下面几个公式来计算：
其中:
最后:
改进算法
FR算法改进了弹簧算法，是现在用途最为广泛的布点算法，很多算法都是在这个算法上改进的。
FR受到了天体重力系统的启发，使用力来计算每个节点的速度，而不是加速度，从而得到每个节点应当移动的距离。它的每次迭代分为三个步骤：
- 计算节点之间的排斥力
- 计算相邻节点之间的吸引力
- 综合吸引力和排斥力，通过最大位移限制移动的距离
使用模拟退火算法，使得在图变得越来越稳定时，温度变得更低，节点每次移动的距离就变得更小。其主要原因是防止震荡。
KK算法使得能量最小化，在图的布局上减少了边的交叉，除了需要计算所有节点对之间的最短路径，并不需要其他理论知识。它虽然每一步的计算复杂度高于FR算法，但迭代次数较少，使其执行速度和效果都比FR好。
实现和性能
上一部分是两个月前写的，今天发现了，想把剩下的补足。
本来想把之前的布点算法的代码抽出来，做成一个独立的工程，做成附件发布出来，奈何之前学弟写的代码实在太乱乱乱了，我改了老半天，依然还有一大堆的问题，渐渐地我失去了耐心。觉得这篇文章要太监了。
说点感性的认识吧。布点算法很重要，因为能直观的看出复杂网络的结构和关系。但在点很多时，性能就得考虑，在时间上，几百个点的计算几乎瞬间即可完成，上万个节点的图，采用高性能的布点算法，在牺牲效果的情况下，能在两分钟内跑完，使用力导引算法，则要花费10分钟以上。这可能和学弟实现的代码不够高效率有关系。
附录给出了一个简单的力导引实现，接口应该很容易，一看就能明白。
目前，推荐SM算法，相关链接：
https://en.wikipedia.org/wiki/Stress_Majorization
不过，不建议重复造轮子，有现成的graphviz可以使用。应该把时间花在更重要的事情上。如果我当时没有对语言间互操作那么反感，就不用耽搁那么多功夫重写C#版本了。
其实原本的代码中，包含了SM,超大规模布点，力导引这三种算法的完整代码，有需要的，可以站内我。如果你看到了这里，说明你还是蛮有耐心的。
附录：C#实现
using System;
using System.Collections.Generic;
using System.Linq;
namespace GraphdotNet.ForceAtlas2Algo
{
internal class ForceAtlas2
{
#region Constructors and Destructors
public ForceAtlas2(int nodeNum)
{
nodes = new Node[nodeNum];
for (var i = 0; i < nodes.Length; i++)
{
nodes[i] = new Node();
}
Degree = new int[nodeNum];
}
#endregion
#region Methods
private double GetEdgeWeight(int edgeIdx)
{
return Adjacency[3*edgeIdx + 2];
}
#endregion
#region Constants and Fields
internal double Speed;
private readonly Node[] nodes;
private double outboundAttCompensation = 1;
private Region rootRegion;
#endregion
#region Properties
internal bool AdjustSizes { get; set; }
internal bool BarnesHutOptimize { get; set; }
internal double BarnesHutTheta { get; set; }
internal double EdgeWeightInfluence { get; set; }
internal double Gravity { get; set; }
internal double JitterTolerance { get; set; }
internal bool LinLogMode { get; set; }
internal bool OutboundAttractionDistribution { get; set; }
internal double ScalingRatio { get; set; }
internal bool StrongGravityMode { get; set; }
private double[] Adjacency { get; set; }
private int[] Degree { get; }
private int[] EdgeAdjacency { get; set; }
#endregion
#region Public Methods
public bool CanAlgo()
{
return (nodes.Length != 0);
}
public void EndAlgo()
{
}
public IEnumerable<double> GetX()
{
return nodes.Select(node => node.X);
}
public IEnumerable<double> GetY()
{
return nodes.Select(node => node.Y);
}
public void GoAlgo()
{
if (nodes.Length == 0)
{
return;
}
////initialize node data
var nodeNum = nodes.Length;
foreach (var n in nodes)
{
n.Olddx = n.Dx;
n.Olddy = n.Dy;
n.Dx = 0;
n.Dy = 0;
}
// If Barnes Hut active, initialize root RegionImpl
if (BarnesHutOptimize)
{
rootRegion = new Region(nodes);
rootRegion.BuildSubRegions();
}
// If outboundAttractionDistribution active, compensate.
if (OutboundAttractionDistribution)
{
outboundAttCompensation = 0;
foreach (var n in nodes)
{
outboundAttCompensation += n.Mass;
}
outboundAttCompensation /= nodeNum;
}
// Repulsion (and gravity)
var repulsion = ForceFactory.Builder.BuildRepulsion(AdjustSizes, ScalingRatio);
const int @from = 0;
var to = nodeNum;
var rgCalculator = new RepulsionGravityCalculate(
nodes,
from,
to,
BarnesHutOptimize,
BarnesHutTheta,
Gravity,
StrongGravityMode ? (ForceFactory.Builder.GetStrongGravity(ScalingRatio)) : (repulsion),
ScalingRatio,
rootRegion,
repulsion);
rgCalculator.Run();
// Attraction
var attraction = ForceFactory.Builder.BuildAttraction(
LinLogMode,
OutboundAttractionDistribution,
AdjustSizes,
1*(OutboundAttractionDistribution ? (outboundAttCompensation) : (1)));
var edgeNum = EdgeAdjacency.Length/2;
const double epsilon = 1e-6;
if (Math.Abs(EdgeWeightInfluence - 0) < epsilon)
{
for (var i = 0; i < edgeNum; i++)
{
attraction.Apply(
nodes[EdgeAdjacency[2*i]], nodes[EdgeAdjacency[2*i + 1]], 1);
}
}
else if (Math.Abs(EdgeWeightInfluence - 1) < epsilon)
{
for (var i = 0; i < edgeNum; i++)
{
attraction.Apply(
nodes[EdgeAdjacency[2*i]],
nodes[EdgeAdjacency[2*i + 1]],
GetEdgeWeight(i));
}
}
else
{
for (var i = 0; i < edgeNum; i++)
{
attraction.Apply(
nodes[EdgeAdjacency[2*i]],
nodes[EdgeAdjacency[2*i + 1]],
Math.Pow(GetEdgeWeight(i), EdgeWeightInfluence));
}
}
// Auto adjust speed
var totalSwinging = 0d; // How much irregular movement
var totalEffectiveTraction = 0d; // Hom much useful movement
foreach (var n in nodes)
{
var swinging = Math.Sqrt((n.Olddx - n.Dx)*(n.Olddx - n.Dx) + (n.Olddy - n.Dy)*(n.Olddy - n.Dy));
totalSwinging += n.Mass*swinging;
// If the node has a burst change of direction, then it's not converging.
totalEffectiveTraction += n.Mass*0.5*
Math.Sqrt(
(n.Olddx + n.Dx)*(n.Olddx + n.Dx) + (n.Olddy + n.Dy)*(n.Olddy + n.Dy));
}
// We want that swingingMovement < tolerance * convergenceMovement
var targetSpeed = JitterTolerance*JitterTolerance*totalEffectiveTraction/totalSwinging;
// But the speed shoudn't rise too much too quickly, since it would make the convergence drop dramatically.
const double maxRise = 0.5; // Start rise: 50%
Speed = Speed + Math.Min(targetSpeed - Speed, maxRise*Speed);
// Apply forces
foreach (var n in nodes)
{
// Adaptive auto-speed: the speed of each node is lowered
// when the node swings.
var swinging = Math.Sqrt((n.Olddx - n.Dx)*(n.Olddx - n.Dx) + (n.Olddy - n.Dy)*(n.Olddy - n.Dy));
//double factor = speed / (1f + Math.sqrt(speed * swinging));
var factor = Speed/(1 + Speed*Math.Sqrt(swinging));
n.X += n.Dx*factor;
n.Y += n.Dy*factor;
}
}
public void InitAlgo(double[] adjacencyInput)
{
foreach (var n in nodes)
{
n.ResetNode();
}
var ran = new Random();
foreach (var n in nodes)
{
n.X = (0.01 + ran.NextDouble())*1000 - 500;
n.Y = (0.01 + ran.NextDouble())*1000 - 500;
}
Speed = 1.0;
Adjacency = adjacencyInput;
var len = Adjacency.Length;
var edgeNum = len/3;
EdgeAdjacency = new int[edgeNum*2];
for (var i = 0; i < edgeNum; i++)
{
EdgeAdjacency[2*i] = (int) Adjacency[3*i];
EdgeAdjacency[2*i + 1] = (int) Adjacency[3*i + 1];
}
for (var i = 0; i < 2*edgeNum; i++)
{
Degree[EdgeAdjacency[i]]++;
}
var nodeNum = nodes.Length;
for (var i = 0; i < nodeNum; i++)
{
nodes[i].Mass = 1 + Degree[i];
}
}
public void OutlierProcess()
{
var maxRadius = double.NegativeInfinity;
foreach (var n in nodes)
{
if (n.Mass >= 2)
{
maxRadius = Math.Max(maxRadius, Math.Sqrt(n.X*n.X + n.Y*n.Y));
}
}
maxRadius += 10;
var maxRadiusCeil = (int) Math.Ceiling(maxRadius);
var ran = new Random();
foreach (var n in nodes)
{
if (n.Mass <= 1)
{
n.X = ran.Next(-maxRadiusCeil*20, maxRadiusCeil*20)/20.0;
n.Y = Math.Pow(-1.0, Math.Floor(n.X))*Math.Sqrt(maxRadius*maxRadius - n.X*n.X);
}
}
}
public void ResetPropertiesValues()
{
var nodeNum = nodes.Length;
// Tuning
ScalingRatio = nodeNum >= 100 ? 2.0 : 10.0;
StrongGravityMode = false;
Gravity = 1.0;
// Behavior
OutboundAttractionDistribution = false;
LinLogMode = false;
EdgeWeightInfluence = 1.0;
// Performance
if (nodeNum >= 50000)
{
JitterTolerance = 10.0;
}
else if (nodeNum >= 5000)
{
JitterTolerance = 1.0;
}
else
{
JitterTolerance = 0.1;
}
BarnesHutOptimize = nodeNum >= 1000;
BarnesHutTheta = 1.2;
}
#endregion
}
}
