在数据可视化的过程中，绘制网络拓扑图是很重要的，它能清晰呈现一个复杂网络的结构，节点的重要性和关系。比如下面几张图：
下面这张图是我的软件绘制的：
这些都有一个共同的问题，就是如何让图绘制的更加美观？ 复杂的图结构，已经没法人工布局了。所以计算机自动布局，就成了一个有趣而重要的话题。我们将其称为布点算法。
基本原理
一个好的布点算法应该能尽量满足以下四个特点：
- 对称性：绘制网络对应将具有相同结构的子图围绕绘图中心平衡布局
- 连接角最大化原则，使同一节点任意两条边形成的角度尽量大
边交叉数量最小原则（最重要）：绘图时应尽量减少相互交叉边的数量。 -
直线边原则：边尽量直线，避免曲线，同时线尽可能短。其中，尽量避免交叉是最重要的，然而，并不是所有的图结构都能满足这个要求。同时，在性能上，还应当满足如下需求：
-
算法应当能处理尽可能复杂和庞大的图，比如上万甚至上百万的节点
- 运算速度还要尽量快，使其能够在短时间内完成操作
- 可交互性，当用户调整了某一个点的位置，或改变了一些参数，算法是否能够动态调节？
如何表示一个网络
图分为有向图和无向图，有权图和无权图。当图有权时，权重反映了节点间的关系，通常值越大，节点关系越紧密。
- 二维矩阵 最简单也最方便的表示方法，当网络结构巨大时，存储空间浪费是惊人的，而且如果图是无向图，则矩阵是对称的。但存取效率最高
- 临接表， 即将边表达为一个数组，
(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 } }