布点算法的原理和实现

---

在数据可视化的过程中，绘制网络拓扑图是很重要的，它能清晰呈现一个复杂网络的结构，节点的重要性和关系。比如下面几张图：

下面这张图是我的软件绘制的：

这些都有一个共同的问题，就是如何让图绘制的更加美观？ 复杂的图结构，已经没法人工布局了。所以计算机自动布局，就成了一个有趣而重要的话题。我们将其称为布点算法。

基本原理

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

• 对称性：绘制网络对应将具有相同结构的子图围绕绘图中心平衡布局
• 连接角最大化原则，使同一节点任意两条边形成的角度尽量大
  边交叉数量最小原则（最重要）：绘图时应尽量减少相互交叉边的数量。

• 直线边原则：边尽量直线，避免曲线，同时线尽可能短。其中，尽量避免交叉是最重要的，然而，并不是所有的图结构都能满足这个要求。同时，在性能上，还应当满足如下需求：

• 算法应当能处理尽可能复杂和庞大的图，比如上万甚至上百万的节点

• 运算速度还要尽量快，使其能够在短时间内完成操作
• 可交互性，当用户调整了某一个点的位置，或改变了一些参数，算法是否能够动态调节？

如何表示一个网络

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

• 二维矩阵 最简单也最方便的表示方法，当网络结构巨大时，存储空间浪费是惊人的，而且如果图是无向图，则矩阵是对称的。但存取效率最高
• 临接表， 即将边表达为一个数组，(x,y,dist) *N 对于稀疏图来说，存储效率很高，但读写效率低。一种简单的优化方法是对x和y 分别进行排序，排序后按二分查找，即可实现快速存取。
• 实时计算 当网络图巨大，而计算关系本身的运算并不复杂时，可以采用，此时dist[x][y] 被转换为一个函数 GetDistance(x,y), 实时完成计算。这是一种很有趣且比较实用的方法，尤其是当网络图稀疏时。
  为了简化，我们采用二维矩阵作为存取结构，相信读者可以很容易的将其改写为其他形式。

算法介绍

 

力导引算法

力导引布局的方法可以产生相当优美的网络布局，并充分展现网络的整体结构及其自同构特征。该方法最早由Eades在1984年提出。
基本思想是将网络看成一个顶点为钢环，边为弹簧的物理系统，系统被赋予某个初始状态以后，弹簧弹力（引力和斥力）的作用会导致钢环移动，这种运动直到系统总能量减少到最小值停止。每次循环都要计算任意一对点的斥力和相邻两个点的引力，计算复杂度为O（N2）。

基本上绝大多数算法都遵循着这样的原则，即：

1. 将网络看成一个顶点为钢环，边为弹簧的物理系统
2. 不断迭代，使整个系统的总能量达到最小

为了简化，弹簧的受力并不遵循胡克定律，而是自己定义的受力公式，同时决定引力只在相邻节点间。因此，最基本的弹簧算法的伪代码如下：

随机分布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
    }
}