布点算法的原理和实现

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

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

布点算法的原理和实现

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

基本原理

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

  • 对称性:绘制网络对应将具有相同结构的子图围绕绘图中心平衡布局
  • 连接角最大化原则,使同一节点任意两条边形成的角度尽量大
    边交叉数量最小原则(最重要):绘图时应尽量减少相互交叉边的数量。
  • 直线边原则:边尽量直线,避免曲线,同时线尽可能短。其中,尽量避免交叉是最重要的,然而,并不是所有的图结构都能满足这个要求。同时,在性能上,还应当满足如下需求:

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

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

如何表示一个网络

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

  • 二维矩阵 最简单也最方便的表示方法,当网络结构巨大时,存储空间浪费是惊人的,而且如果图是无向图,则矩阵是对称的。但存取效率最高
  • 临接表, 即将边表达为一个数组,(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
    }
}

上一篇:Python3在指定路径下递归定位文件中出现的字符串


下一篇:开发者社区精选直播合集(十七)| 云计算十年的发展与应用