基本原理

一个好的布点算法应该能尽量满足以下四个特点：

对称性：绘制网络对应将具有相同结构的子图围绕绘图中心平衡布局
连接角最大化原则，使同一节点任意两条边形成的角度尽量大
边交叉数量最小原则（最重要）：绘图时应尽量减少相互交叉边的数量。
直线边原则：边尽量直线，避免曲线，同时线尽可能短。其中，尽量避免交叉是最重要的，然而，并不是所有的图结构都能满足这个要求。同时，在性能上，还应当满足如下需求：
算法应当能处理尽可能复杂和庞大的图，比如上万甚至上百万的节点
运算速度还要尽量快，使其能够在短时间内完成操作
可交互性，当用户调整了某一个点的位置，或改变了一些参数，算法是否能够动态调节？

如何表示一个网络

图分为有向图和无向图，有权图和无权图。当图有权时，权重反映了节点间的关系，通常值越大，节点关系越紧密。

二维矩阵 最简单也最方便的表示方法，当网络结构巨大时，存储空间浪费是惊人的，而且如果图是无向图，则矩阵是对称的。但存取效率最高
临接表，即将边表达为一个数组，(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 Destructorspublic 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 Methodsprivate double GetEdgeWeight(int edgeIdx){return Adjacency[3*edgeIdx + 2];}#endregion#region Constants and Fieldsinternal double Speed;private readonly Node[] nodes;private double outboundAttCompensation = 1;private Region rootRegion;#endregion#region Propertiesinternal 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 Methodspublic 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 datavar 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 RegionImplif (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();// Attractionvar 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 speedvar totalSwinging = 0d; // How much irregular movementvar totalEffectiveTraction = 0d; // Hom much useful movementforeach (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 * convergenceMovementvar 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 forcesforeach (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;// TuningScalingRatio = nodeNum >= 100 ? 2.0 : 10.0;StrongGravityMode = false;Gravity = 1.0;// BehaviorOutboundAttractionDistribution = false;LinLogMode = false;EdgeWeightInfluence = 1.0;// Performanceif (nodeNum >= 50000){JitterTolerance = 10.0;}else if (nodeNum >= 5000){JitterTolerance = 1.0;}else{JitterTolerance = 0.1;}BarnesHutOptimize = nodeNum >= 1000;BarnesHutTheta = 1.2;}#endregion}
}