布點算法的原理和實現

---

在數據可視化的過程中，繪制網絡拓撲圖是很重要的，它能清晰呈現一個復雜網絡的結構，節點的重要性和關系。比如下面幾張圖：

下面這張圖是我的軟件繪制的：

這些都有一個共同的問題，就是如何讓圖繪制的更加美觀？ 復雜的圖結構，已經沒法人工布局了。所以計算機自動布局，就成了一個有趣而重要的話題。我們將其稱為布點算法。

基本原理

一個好的布點算法應該能盡量滿足以下四個特點：

• 對稱性：繪制網絡對應將具有相同結構的子圖圍繞繪圖中心平衡布局
• 連接角最大化原則，使同一節點任意兩條邊形成的角度盡量大
  邊交叉數量最小原則（最重要）：繪圖時應盡量減少相互交叉邊的數量。

• 直線邊原則：邊盡量直線，避免曲線，同時線盡可能短。其中，盡量避免交叉是最重要的，然而，並不是所有的圖結構都能滿足這個要求。同時，在性能上，還應當滿足如下需求：

• 算法應當能處理盡可能復雜和龐大的圖，比如上萬甚至上百萬的節點

• 運算速度還要盡量快，使其能夠在短時間內完成操作
• 可交互性，當用戶調整了某一個點的位置，或改變了一些參數，算法是否能夠動態調節？

如何表示一個網絡

圖分為有向圖和無向圖，有權圖和無權圖。當圖有權時，權重反映了節點間的關系，通常值越大，節點關系越緊密。

• 二維矩陣 最簡單也最方便的表示方法，當網絡結構巨大時，存儲空間浪費是驚人的，而且如果圖是無向圖，則矩陣是對稱的。但存取效率最高
• 臨接表， 即將邊表達為一個數組，(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
    }
}