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)
ノード 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) # ノード番号をつけて、ノードの配置は円形
ノードの数とエッジの数だけ指定して、つながり方はランダムに決めたい場合は次のようにすればOK。
numnodes, numedges = 10, 20 G = Graph(numnodes, numedges)
有向グラフにしたい場合は、Graph
をDiGraph
にすればよい。
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)
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)
分布が指数関数的に減衰しているので良さげ。