Directed Scale-Free Graphsによるグラフの生成 - NO!と言えるようになりたい

Webリンクなどは有向グラフとよばれるグラフとなっており，リンクに向きがあります．Web用のソフトウェアを作成したとき，シミュレーション等を行うために疑似データを生成したくなりますが，"Directed scale-free graphs", B. Bollobas, C. Borgs and O. Riordan, Proceedings of the 14th Annual ACM-SIAM Symposium on Discrete Algorithms (SODA), 132-139 (2003) (http://research.microsoft.com/~chayes/)に示されている方式を用いると，いわゆるWebのようなつながりを持つグラフを生成する事が出来ます．これは名前の通り，スケールフリー性をもっているため，いろいろな用途に使えそうです．

アルゴリズムは少し複雑です．

アルゴリズムの説明に移る前に，二つの確率密度関数を定義しておきます．

Pin(v = vi) = (din(vi) + δin) / (t + δin・n(t))

と

Pout(v = vi) = ((dout(vi) + δout) / (t + δout・n(t))

です．

ただしここで，din(vi) と dout(vi) はそれぞれノード vi に入ってくるリンクの数と出て行くリンクの数となり，t は時刻（リンクの合計数でもある），n(t) は時刻 t でのノードの合計となります．

δin と δout はどのようにリンクを張るかのパラメータとなり，定数となります．関数 Pin() は入ってくるリンクが多いほど，リンクを張られる確率が高くなっていますが，δin を大きくすると，その確率の差を小さくする事になります．δout も同じ働きです．

あとは，下記に記されたステップを繰り返すと生成する事が出来ます．

確率αで新たなノード v を生成し，その v から既存のノード w へとリンクを張る．ただし，w は確率密度関数 Pin() に従ってランダムに選ばれる．
確率βで既存のノード v と w を選び，v から w へとリンクを張る．ただし，v, w はそれぞれ，確率密度関数 Pout(), Pin() に従ってランダムに選ばれる．
確率 1 - α - βで新たなノード w を生成し，既存のノード v から w へとリンクを張る．ただし，v は確率密度関数 Pout() に従ってランダムに選ばれる．

Bollobasらの論文によると，Webと同じような特徴を持つグラフを生成したい場合，α = 0.41，β = 0.49，δin = 0，δout = 0，もしくは，α = 0.41，β = 0.59，δin = 0.24，δout = 0 とすればよいらしいです．

Pythonでのプログラムは以下の通りとなります．

#!/usr/bin/env python
# this program generates directed scale-free graphs
#
# see "Directed scale-free graphs", B. Bollobas, C. Borgs and O. Riordan,
# Proceedings of the 14th Annual ACM-SIAM Symposium on Discrete Algorithms
# (SODA), 132-139 (2003)

import random


class node:
    def __init__(self, id):
        self.inLink  = set()
        self.outLink = set()
        self.id = id


class dirsf:
    def _initialize(self):
        random.seed()

        self.nodes = {}
        self.nodes[0] = node(0)
        self.nodes[1] = node(1)
        self.nodes[2] = node(2)

        self.nodes[0].outLink.add(1)
        self.nodes[1].inLink.add(0)
        self.nodes[1].outLink.add(2)
        self.nodes[2].inLink.add(1)
        self.nodes[2].outLink.add(0)
        self.nodes[0].inLink.add(2)

        self.edge = 3
        self.n    = 2


    def generate(self, step, alpha, beta, din, dout):
        self._initialize()

        for var in range(step):
            rnd = random.random()
            print self.n

            if rnd < alpha:
                self.n += 1
                fromNode = node(self.n)
                toNode   = self._getNode(din, self._probIn)

                fromNode.outLink.add(toNode.id)
                toNode.inLink.add(fromNode.id)

                self.nodes[self.n] = fromNode
            elif rnd < alpha + beta:
                while True:
                    fromNode = self._getNode(dout, self._probOut)
                    toNode   = self._getNode(din,  self._probIn)

                    if toNode.id in fromNode.outLink:
                        continue
                    
                    if fromNode.id != toNode.id:
                        fromNode.outLink.add(toNode.id)
                        toNode.inLink.add(fromNode.id)
                        break
            else:
                self.n += 1
                fromNode   = self._getNode(dout, self._probOut)
                toNode = node(self.n)

                fromNode.outLink.add(toNode.id)
                toNode.inLink.add(fromNode.id)

                self.nodes[self.n] = toNode

            self.edge += 1


    def _probIn(self, n, delta):
        din = float(len(self.nodes[n].inLink))
        return (din + delta) / (self.edge + delta * len(self.nodes))


    def _probOut(self, n, delta):
        dout = len(self.nodes[n].outLink)
        return (dout + delta) / (self.edge + delta * len(self.nodes))


    def _getNode(self, delta, prob):
        while True:
            n = int(random.uniform(0, len(self.nodes)))
            rnd = random.random()
            p = prob(n, delta) / len(self.nodes)

            if rnd < p:
                return self.nodes[n]


    def printLink(self):
        for n in range(self.n):
            for node in self.nodes[n].outLink:
                print '(%d, %d)' % (n, node)


    # draw with NodeBox
    # see http://nodebox.net/code/index.php/Graph
    def draw(self, width, height):
        try:
            graph = ximport("graph")
        except ImportError:
            graph = ximport("__init__")
            reload(graph)

        size(width, height)
        g = graph.create()

        for n in range(self.n):
            if len(self.nodes[n].outLink) == 0:
                continue

            for id in self.nodes[n].outLink:
                g.add_edge(n, id)

        g.solve()
        g.styles.apply()
        g.draw(directed = True, traffic = 1)
            

if __name__ == "__main__":
    graph = dirsf()
    graph.generate(1000, 0.41, 0.59, 0.24, 0.0)

    graph.printLink()

"""
graph = dirsf()
graph.generate(1000, 0.41, 0.59, 0.24, 0.0)

graph.draw(1440, 900)
"""

NodeBoxを用いてvisualizeしてみたら，以下のようになりました．