Julia - ネットワーク上(グラフ上)のランダムウォーク シミュレーション

f:id:goropikarikun:20171028012853g:plain:w300

Julialangを使ってグラフ上のランダムウォークのシミュレーションをしてみる。
今回、グラフにはLightGraphs.jlを使用した。

環境

julia> versioninfo()
Julia Version 0.6.0
Commit 903644385b* (2017-06-19 13:05 UTC)
Platform Info:
  OS: Linux (x86_64-pc-linux-gnu)
  CPU: Intel(R) Core(TM) i5-4460T CPU @ 1.90GHz
  WORD_SIZE: 64
  BLAS: libopenblas (HASWELL)
  LAPACK: liblapack
  LIBM: libm
  LLVM: libLLVM-3.9.1 (ORCJIT, haswell)


julia> println("LightGraphs.jl v", Pkg.installed("LightGraphs"))
LightGraphs.jl v0.11.0

julia> println("GraphPlot.jl v", Pkg.installed("GraphPlot"))
GraphPlot.jl v0.2.0

サンプル

ものすごくシンプルな場合だったら以下のようなコードで終わる。

julia> nnodes, nedges = 10, 20
julia> G = Graph(nnodes, nedges)
julia> gplot(G, nodelabel=1:nv(G), layout=circular_layout)

f:id:goropikarikun:20171028015735p:plain:w300

ノード 1 から 10 ステップ分歩かせる。

julia> startnode, niter = 1, 10; randomwalk(G, startnode, niter)
10-element Array{Int64,1}:
  1
  4
  1
  5
  7
 10
  9
 10
  8
  4

パッケージをインストール

Pkg.add("LightGraphs")
Pkg.add("GraphPlot")
using LightGraphs, GraphPlot

グラフを作る

ランダムウォーカーが動き回るためのグラフをまず作る。

G = Graph() # 無向グラフ
add_vertices!(G, 10)  # ノードを10個作る
add_edge!(G, 1, 10) # ノード1とノード10をつなぐ
add_edge!(G, 2, 8) # ノード2とノード8をつなぐ
# ... 以下同様

gplotで作ったグラフを図示することができる。

gplot(G, nodelabel=1:nv(G), layout=circular_layout) # ノード番号をつけて、ノードの配置は円形

f:id:goropikarikun:20171028022805p:plain:w300

ノードの数とエッジの数だけ指定して、つながり方はランダムに決めたい場合は次のようにすればOK。

numnodes, numedges = 10, 20
G = Graph(numnodes, numedges)

有向グラフにしたい場合は、GraphDiGraphにすればよい。

Random walk

LightGraphs.jlには最初からランダムウォークのための関数であるrandomwalkという名前そのままの関数がある。
始点とステップ数を指定するとランダムウォーカーの軌跡が返ってくる。

snode, niter = 1, 10
randomwalk(G, snode, niter)

10-element Array{Int64,1}:
 1
 3
 4
 1
 3
 1
 3
 4
 2
 1

初到達時間 - First passage time(FPT)

実際に研究で使おうと思うとFPTが知りたくなるので、元のランダムウォークのコード*1 を改変しFPTを調べる関数を定義する。

"""
    fpt(g, s, f ;niter)
First passsage time

Perform a random walk on graph `g` starting at vertex `s` and finishing at vertex `f`.
niter is upper bound of the number of iterations to avoid infinite loop. Default is 10*nv(g)
Return a vector of vertices visited in order and total time step.
"""
function fpt(g::AG, s::Integer, f::Integer ;niter::Integer=10*nv(g)) where AG <: AbstractGraph{T} where T
    s in vertices(g) || throw(BoundsError())
    f in vertices(g) || throw(BoundsError())
    Bool(in.(s, connected_components(g))' * in.(f, connected_components(g)))|| error("Node $s and node $f are disconnected")
    visited = Vector{T}()
    sizehint!(visited, niter)
    currs = s
    i = 1
    while true
        push!(visited, currs)
        i += 1
        currs = rand(out_neighbors(g, currs))
        if currs == f
            push!(visited, currs)
            break
        end
    end
    return visited[1:i], i-1
end
V, E = 10, 20
G =  Graph(V, E)
gplot(G, nodelabel=1:V, layout=circular_layout)

f:id:goropikarikun:20171028103630p:plain:w300

srand(2017); 
fptresult = fpt(G, 1, 1) # ノード1からスタートしてノード1に戻ってくるまでの時間

([1, 4, 3, 4, 9, 5, 1], 6)

あるノードからスタートしてもう一度スタート地点に戻ってくるのにかかる時間の確率分布$$p_v(t)$$をプロットしてみる。

using Plots; gr()
@time begin
    srand(1)
    g = Graph(200, 10000)
    node = 4
    @show mean(degree(g))
    @show degree(g, node)
    nsample = 1_000
    first_passage_time = ones(Int, nsample)
    for i in 1:nsample
        first_passage_time[i] = fpt(g, node, node)[2] # first passege time of 'node' to 'node'
    end
end
histogram(first_passage_time, yscale=:log10, norm=true, leg=false, xlabel="time step", ylabel="probability distribution")


mean(degree(g)) = 100.0
degree(g, node) = 108
  0.121123 seconds (94.16 k allocations: 50.032 MiB, 10.71% gc time)

f:id:goropikarikun:20171028105429p:plain

分布が指数関数的に減衰しているので良さげ。