MATLAB移民のためのJulia tips

MATLAB移民のためのJulia tips

何も考えずに速く計算できないのならば、何もやりたくない。

Indicator行列：broadcast

サンプル

Indicator行列

例えば、各要素は１〜３番目のどれかが１、他は０になっていて、 要素１：１番目が１、他は０
要素２：３番目が１、他は０
要素３：２番目が１、他は０
要素４：１番目が１、他は０
というとき、
1 0 0
0 0 1
0 1 0
1 0 0
という行列。

作り方：broadcast()する

上記のindicator行列を作るとき、例えば配列として

v = [1, 3, 2, 1] #ベクトルの場合はvec(v)

を持っていたとすると、次のようにして作れる*1：

IND = broadcast(.==, v, (1:3)')

4x3 Array{Int64,2}:
1 0 0
0 0 1
0 1 0
1 0 0

  • (1:3)'はベクトル[1 2 3]。
  • vと(1:3)'を要素ごとに評価（.==）し、値が等しいところで１が立つ。

Boradcastについては、以下参照：
Multi-dimensional Arrays — Julia Language 0.4.3-pre documentation （本家） NumPyのブロードキャストのメモ - 唯物是真 @Scaled_Wurm（Pythonの記事。わかりやすいです。）

*1:ここでの内容はstackoverflow julia lang - convert Array to indicator matrix - Stack Overflow のajcr氏の回答の内容ほぼそのままです。

プロット：Gadflyのtips

プロット

Juliaの上でプロットしたかったら、PyPlotかGadflyが標準的らしい。

  • Gadflyは配色のセンスがいい。
  • 他のに比べて重い。（特に最初に読み込むときは信じがたいほど遅い）

もっと比べたい人はこちらへ: Graphics in Julia

DataFrameを使ってプロットする方法と、配列や関数をそのままプロットする方法の２通りあるらしい。

  • 素朴なプロットならDataFrameなしの方がお手軽。
  • ちょっと複雑なプロットならDataFrameを使った方がよい。

基本DataFrame推しのようなので、諦めてDataFrameを使った書き方に慣れた方がよさそう。

準備

JuliaのターミナルでGadfly, DataFramesを入れる。

Pkg.add("Gadfly")
Pkg.add("DataFrames") # <-複数形ね

(Cairoも入れておくとよいと言われているけど、今のところCairoで解決した問題がないので保留 [1/2/2016])
使うときは、最初に

using Gadfly
using DataFrames # <-複数形ね

と書いてから始める。 (欲張って using PyPlotも書くと、競合して動かない(っぽい)ので注意。)

配列データをプロットする

２次元density plot: rectbin

x-座標(x)、y-座標(y)、(x,y)での値(value)の３つを、ベクトルX, Y, Zとして与える：

plot(x=X,y=Y,color=Z,Geom.rectbin)

<例> N-by-M行列PSIの行列要素をdensity plotする。

X = Int64[]
Y = Int64[]
Z = Float32[]
for y = 1:size(PSI,1)
    for x = 1:size(PSI,2)
        push!(X,x)
        push!(Y,y)
        push!(Z,PSI[y,x])
    end
end

plot(x=X,y=Y,value=Z,Geom.rectbin)

f:id:tatmann9:20160102213832p:plain

ref: (参考にさせていただいた記事が見つからない。。。)

DataFrameを使ってプロットする

こんなサンプルデータを考える*1。

n  = 10
x  = collect(1:n)
y1 = rand(n)
y2 = rand(n)
y3 = rand(n)

配列データをDataFrameにする

# Put the data in a DataFrame
using DataFrames
d = DataFrame( 
  x = vcat(x,x,x),
  y = vcat(y1,y2,y3),
  group = vcat( rep("1",n), rep("2",n), rep("3",n) )
)
  • データは行列ではなく、配列にしておかないとダメ。行列になっていたらvec()で変換。
  • repはreplicateの略で、n成分の配列を作ってくれる*2。
  • x, y, groupは名前を変えてもよい。これで軸名を指定できる。

scatter plot

# Plot
using Gadfly
plot( 
  d, 
  x=:x, y=:y, color=:group, 
  Geom.point,
  Scale.discrete_color_manual("green","red","blue")
)
  • Scale.discrete_color_manual("green","red","blue")はいろいろwarningが出るけど、とりあえず実行できた。
  • colorは名前を変えてはダメ。ぜったい。

histogram

単純。

using DataFrames
using Gadfly
v = [1,2,2,0,0,0,4,3]  # 例えばこんなデータを考える
d = DataFrame( 
    x = v # vが配列データじゃなかったら、vec(v)で配列にする
)
plot(d, x="x", Geom.histogram(bincount=6))
  • オプションのbincountはなくてもOK。
  • 整数要素の場合、x=0のbinがバグるっぽいので、微小乱数でも付けておけば大丈夫かと。（整数だと、常にbinの境界の値を判定しているから。たぶん。）
  • 「規格化のオプションがない！」と思ったら、それはGeom.histogramではなく、Geom.densityを使えばよい。

プロットに軸名、タイトルをつける

例としてrectbinのところで用いた行列のdensity plotを扱う。
x軸とy軸の軸名、プロットのタイトルをつけるには、

plot(
x=X,y=Y,color=Z,
Guide.xlabel("row"), 
Guide.ylabel("column"), 
Guide.title("matrix"),
Geom.rectbin)

でできる*3。

随時更新予定

(以上、R移民にしてみれば常識かもしれないが、MATLAB移民に言わせればry)

*1: 以下の数小節は、julia lang - How to overplot multiple arrays of data in one Gadfly plot? - Stack Overflow のVincent Zoonekynd氏の回答の内容ほぼそのままです。

*2:Rプログラム (TAKENAKA's Web Page)

*3:Introducing Julia/Plotting - Wikibooks, open books for an open world

型：配列と行列

基本操作

「一次元配列も1-by-N行列も似たようものだから、どっちでも動くはず」と思ってハマっているMATLAB移民はこちら。
普通の言語は、そのへん適当にできるようにはなっていない。

1-by-N行列を一次元配列に変換

vec()で1-by-N行列を一次元配列にできる。

u = [1 2 3] # 1x3 Array{Int64,2}
vec(u) # 3-element Array{Int64,1}

一次元配列を1-by-N行列に変換

転置操作'で一次元配列を1-by-N行列にできる。
２回転置すればN-by-1行列にできる。

u = collect(1:3) # 3-element Array{Int64,1}
v = u' # 1x3 Array{Int64,2}
v = u'' # 3x1 Array{Int64,2}

（もっとスマートな書き方あれば教えてください。）

一応書いておくと

a = [1 2 3]とすると、1-by-3の行列となり、
a = [1, 2, 3]は、一次元配列。
ones(3)とすると、[1.0, 1.0, 1.0]の配列になるが、
ones(1,3)は[1.0 1.0 1.0]の1-by-3の行列となる。

深さ優先探索(DFS)：グラフの連結成分

サンプル

どこにでも書いてある、DFSで連結成分を出してくるという操作のサンプル。
特にJuliaだから特殊ということはないはず。
枝の集合linksが与えられたときに、連結成分の頂点CCnodesを返す。

Stackでやる方法

function dfs_stack(nb,root)
    visited = Int64[]
    stack = push!(Int64[],root)
    while !isempty(stack)
        node = pop!(stack)
        if node in visited
            continue
        else
            push!(visited,node)
            append!(stack,filter(x->!(x in visited), nb[node]))
        end
    end
    return visited
end

links = [1 2; 2 3; 2 4; 4 1; 1 5; 4 5;3 4; 7 8; 8 9; 9 7; 4 8]
Ntot = maximum(links)
links = vcat(hcat(links[:,2],links[:,1]),links)
nb = [links[links[:,1].==lnode,2] for lnode = 1:Ntot]

CCnodes = dfs_stack(nb,1)
sort(CCnodes)

8-element Array{Int64,1}:
1
2
3
4
5
7
8
9

再帰でやる方法

function DFS(nb, root, visited=nothing) 
    if visited == nothing
        visited = Int64[]
    end

    if root in visited
        return
    end
    
    push!(visited,root)
    
    for vertex in filter(x->!(x in visited), nb[root])
        DFS(nb, vertex, visited)
    end
    return visited
end

links = [1 2; 2 3; 2 4; 4 1; 1 5; 4 5;3 4; 7 8; 8 9; 9 7; 4 8]
Ntot = maximum(links)
links = vcat(hcat(links[:,2],links[:,1]),links)
nb = [links[links[:,1].==lnode,2] for lnode = 1:Ntot]

CCnodes = DFS(nb,1,nothing)
sort(CCnodes)

結果は前と同じ。

引数は参照渡し（代入で困った人はこちら）

基本操作

例えば、こういうこと。

a = [1,2,3]
b = a
b[1] = 10
println(b) # => [10,2,3]
println(a) # => [10,2,3]

aをbに代入してから、bの要素を書き変えると、aまで変更されちゃうということ。
参照渡しというやつ。（詳しいことはpythonで参照渡しの説明をしている記事を読んでください。）

ちなみに引数でいじらなければ、大丈夫。

a = [1,2,3]
b = a
b = [10,2,3]
println(b) # => [10,2,3]
println(a) # => [1,2,3]

最初の例も、b = a[:]としておけば、

a = [1,2,3]
b = a[:]
b[1] = 10
println(b) # => [10,2,3]
println(a) # => [1,2,3]