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)
add_edge!(G, 1, 10)
add_edge!(G, 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, 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]
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)
分布が指数関数的に減衰しているので良さげ。