新しいブログを作りましたのでお知らせします
新しいブログはhttps://www.bicycle1885.org/です。 こちらは今後更新しない予定です。 よかったら最初の記事をご覧ください。 www.bicycle1885.org
Juliaのfor文はなぜ速い?
1年以上ぶりのブログ更新ですね。
atsuoishimotoさんのブログ記事で,私のブログ記事が言及されていました。 atsuoishimoto.hatenablog.com
書いてある内容はまさにその通りで,for文を回すこと自体が取り立てて遅いわけではなく,その中で何らかの処理をする際に動的型付け言語であるPythonの性質上避けて通れない処理があり,その処理をfor文の中で何度も繰り返すことで結果としてfor文全体の処理が遅くなるという内容でした。
そうすると,ひとつの疑問として,なぜ同じ動的型付け言語であるはずのJulia言語では同様の処理が高速なのかが気になるところでしょう。この記事ではそのあたりを少し解説しようと思います。
Juliaの型推論
問題となっているコードは,以下のように2つのベクトルの内積を計算するコードです。
function dot(a, b) s = zero(eltype(a)) for i in 1:endof(a) s += a[i] * b[i] end return s end
このコードは,特に型の指定などがないにも関わらず,Pythonで似たコードを書いた場合と比べて非常に高速に動作します。
Julia:
julia> a = ones(100000); b = ones(100000); julia> @benchmark dot(a, b) BenchmarkTools.Trial: memory estimate: 16 bytes allocs estimate: 1 -------------- minimum time: 121.605 μs (0.00% GC) median time: 121.666 μs (0.00% GC) mean time: 132.996 μs (0.00% GC) maximum time: 1.456 ms (0.00% GC) -------------- samples: 10000 evals/sample: 1
In [2]: %timeit dot(a, b) 31.4 ms ± 680 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)
ここでPythonがこれほど遅くなる理由は,元のブログ記事で解説されているとおり,for文の中で行われている一つ一つの計算の裏で様々な付随する処理が行われているためです。例えば,乗算を行う*
演算子では,そもそも被演算子の型が乗算をサポートしているかなどの確認が必要です。動的言語では変数の型が動的に決定されるため,こうした処理が付随するのは避けて通れないようにも思えますが,Juliaはこうした処理を省略することで高速化を達成しています。
Juliaはなぜこうした処理を省略できるのでしょうか? その答えはJuliaの型推論にあります。実は,Juliaの処理系は関数を実行する前に型推論を行います。Juliaの型推論器は,関数のコードを見て回って,型が決定できるところに自動的に型注釈を付けて回ります。先ほどのコードを例に少し説明してみましょう。このコード例では,dot(a, b)
という呼び出しが起きた際に引数であるa
とb
の型が決定します。上の例では両方ともVector{Float64}
(これはArray{Float64,1}
の別名) というFloat64
型 (倍精度の浮動小数点数の型) を要素とする1次元配列です。こうすると当然for文の中にあるa[i]
やb[i]
はFloat64
型の値を取り出すことが分かります。すると,a[i] * b[i]
という配列のi
番目の要素を取り出して積を計算する処理は,Float64
同士の積だと分かりますので,被演算子が乗算をサポートしているかどうかをいちいちfor文の内部で確認する必要はありません。こうした処理を勝手に行ってくれるのが型推論です。
型推論の結果を見てみましょう。@code_typed
というマクロを使うと,関数の型推論の結果を表示することができます。
julia> @code_typed dot(a, b) CodeInfo(:(begin s = (Base.sitofp)(Float64, 0)::Float64 # line 3: $(Expr(:inbounds, false)) # 省略 SSAValue(4) = (Base.arraysize)(a, 1)::Int64 # 省略 $(Expr(:inbounds, :pop)) SSAValue(5) = (Base.select_value)((Base.slt_int)(SSAValue(4), 0)::Bool, 0, SSAValue(4))::Int64 SSAValue(6) = (Base.select_value)((Base.sle_int)(1, SSAValue(5))::Bool, SSAValue(5), (Base.sub_int)(1, 1)::Int64)::Int64 #temp# = 1 17: unless (Base.not_int)((#temp# === (Base.add_int)(SSAValue(6), 1)::Int64)::Bool)::Bool goto 27 SSAValue(7) = #temp# SSAValue(8) = (Base.add_int)(#temp#, 1)::Int64 i = SSAValue(7) #temp# = SSAValue(8) # line 4: s = (Base.add_float)(s, (Base.mul_float)((Base.arrayref)(a, i)::Float64, (Base.arrayref)(b, i)::Float64)::Float64)::Float64 25: goto 17 27: # line 6: return s end))=>Float64
下の方を見てみると,(Base.mul_float)((Base.arrayref)(a, i)::Float64, (Base.arrayref)(b, i)::Float64)::Float64
と書かれた部分があるのがわかると思います。この部分がまさにa[i] * b[i]
の計算にに対応する部分で,それぞれの被演算子と演算結果がFloat64
になることがちゃんと推論できています。また,一番最後に=>Float64
と書かれているように,dot(a, b)
の計算結果がFloat64
型であることも推論されています。計算結果の型も推論できることで,他の関数からdot
関数を呼び出したときもその計算結果が推論でき,型の情報がどんどん伝播していきます。
型推論は,ソースコードの上でインタープリタを実行するように行います。このインタープリタは引数の値などを扱うインタープリタではなく,型を扱うインタープリタです。型推論の詳しいアルゴリズムは実装を読むか,Inference Convergence Algorithm in Juliaと Inference Convergence Algorithm in Julia - Revisitedなどの解説を参照して下さい。
コードの生成
ソースコードに現れる変数や関数の返り値の型が型推論で判明すると,いくつかの確認処理を省くだけでなく,Juliaはその型に特化した命令を生成することができます。関数の実行時にはコード生成を行うことで処理速度を上げています。その中間状態を覗く@code_llvm dot(a, b)
というマクロがありますので,結果を見てみましょう。
julia> @code_llvm dot(a, b) define double @julia_dot_62600(i8** dereferenceable(40), i8** dereferenceable(40)) #0 !dbg !5 { # 省略 idxend2: ; preds = %idxend %18 = add i64 %"#temp#.08", 1 %19 = getelementptr double, double* %10, i64 %13 %20 = load double, double* %19, align 8 %21 = getelementptr double, double* %12, i64 %13 %22 = load double, double* %21, align 8 %23 = fmul double %20, %22 %24 = fadd double %s.09, %23 %25 = icmp eq i64 %"#temp#.08", %4 br i1 %25, label %L27.loopexit, label %if }
この結果はLLVM IRという言語で書かれています。LLVM IRは機械語に近い低レベルな命令列で,C言語などのコンパイラであるclangやRustのバックエンドでも使われているものです。
この結果を見ると%20 = load double, double* %19, align 8
という命令でメモリのある部分からdouble
の値を読み出しているのがわかります。これはちょうどa[i]
(もしくはb[i]
)に相当する処理です。また,読み出した値の積を計算する%23 = fmul double %20, %22
という命令があるのも分かります。fmul
という命令は浮動小数点数の積を計算するのに特化した命令です ('fmul' Instruction)。ではa
やb
の要素がFloat64
でなく整数型のInt64
だったらどうなるでしょうか? そういう場合でもJuliaは型推論を行ってInt64
型に特化した積の演算命令を生成します。このように型推論の結果を使ってある処理に特化した命令を使うことで,非常に高効率なコードを生成します。ループ1回あたりでは大した違いはありませんが,長く回すループだと積もり積もって大きな差になります。こうして生成された関数のコードは,次回以降に同じ型の引数がくると再利用されます。
副作用
Juliaは関数の実行時に型推論とコード生成を行うことを説明しましたが,これには副作用もあります。初めて関数を呼び出す際に型推論とコード生成・最適化にかかるコストが上乗せされます。このコストがペイするほどのものかどうかを確認すべきなのでしょうが,現在のJulia (v0.6.2)は特殊なことをしない限り必ず型推論とコード生成を行います。ですので,初回の関数呼び出しでワンテンポ遅れたり,それが巨大なパッケージになると積もり積もって数秒から数十秒の遅れになったりします。とても良い例えをTwitterで見つけたので引用しようと思います。
これは人間と車が一メートル競走でどっちが早いかを見てるのと同じことだとおもいます.
— ごまふあざらし (@MathSorcerer) 2018年1月7日
将来的にもっとJuliaの処理系が賢くなれば,自分が今1メートル走をしようとしてるのかフルマラソンをしようとしてるのか分かるようになり,こうした副作用の効果も軽減されるようになると思います。
Pythonで同じことができないのか
Numbaというプロジェクトは,Juliaに似たような型推論とコード生成を行っているようです。
使い方は簡単で,@jit
というアノテーションを関数につけるだけです。
import numba @numba.jit def dot(a, b): s = 0.0 for i in range(len(a)): s += a[i] * b[i] return s
実際にこれを試してみると,Juliaと同じレベルの高速なコードになりました。
In [2]: %timeit dot(a, b) 133 µs ± 1.43 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)
どの程度複雑なコードまで上手く行くのかは分かりません(例えば,ユーザ定義型でもちゃんと最適化できるのか)が,非常に面白いプロジェクトだと思います。今後,RやRubyでも似たようなアプローチが出てくるかもしれませんね。
EzXML.jlを作った話
先月ですが、EzXML.jlというXMLやHTMLを扱うためのパッケージをリリースしました。 EzXML.jlはlibxml2というC言語で書かれたメジャーなライブラリのラッパーなのですが、それをJuliaでラップする過程で色々と工夫がありましたので、その辺を共有したいと思います。また、パッケージの開発についても私の思うところを色々と紹介していきたいとも居ますので、これからJuliaのパッケージを作ってみようと思う方には参考になるかと思います。
何故作ったか
そもそも、JuliaにはLightXML.jlというパッケージがあり、こちらもlibxml2をラップしたパッケージでEzXML.jlと重複する部分がかなりあります。LightXML.jlはJuliaのXMLを扱うパッケージの中では最もよく使われているもので、多くの依存パッケージがあります。しかしながら、いくつか問題があり、私はこのパッケージにパッチを送るでもなく新しいパッケージを作ることにしました。
LightXML.jlで問題だと思われたものを列挙すると主に以下の様なものが挙げられます。 * APIがJuliaっぽくない * 型の設計があまりうまくない * 機能が少ない * C側で確保されたメモリーを自動的に解放してくれない
まず、APIに関してですが妙な名前の関数が提供されています。例えば、XMLファイルをパースして読み込む関数はparse_file
という何をするのか分からない名前です。Juliaではしょっちゅうusing
を使ってパッケージから提供されている関数をパッケージ名を付けずに利用しますから、例えばちょっと大きめのプログラムでparse_file
という関数が呼ばれていても、これがXMLを読み込む関数だとはなかなか思えません。また、幾つかの関数名にはchild_nodes
やis_elementnode
などアンダースコアが使われていて、これもJuliaっぽい関数の命名法ではないです。例えば、Juliaの標準関数を見てみるとeachline
やhaskey
などアンダースコアを使わない関数名が一般的で、スタイルガイドにもそのように記述されています。
型に関しても、ちょっと設計がよくない気がします。XMLのDOMモデルでは、XMLの文書中を構成する要素をノードと呼び、<element/>
のような要素ノードや<!-- comment -->
のようなコメントノードなど様々な種類のノードがあるのですが、LightXML.jlではそのようなノードを別々のJuliaの型でラップしています。ですので、以下のようにノードを見つけるとその種類を判別して目的の型のノードのキャストする必要があります。
for c in child_nodes(xroot) if is_elementnode(c) e = XMLElement(c) # ... end end
これはいかにも面倒くさいですし、一度生成したオブジェクトを別のオブジェクトへキャストすることで余計なオブジェクトの生成も起きますから、パフォーマンス上もあまりよろしくないと思います。
また、libxml2は名前空間やXPathなど様々な機能を実装しているのですが、LightXML.jlから提供されているのはそのごく一部でXMLの読み書きとDOM操作だけです。それに、libxml2を通して確保したメモリー領域は、Julia側からの参照がなくなっても自動的には解放してくれないので、自分で解放しないとメモリーリークを起こします。
これらの問題を解決して快適にXMLを扱えるようにするため、EzXML.jlという新しいパッケージを開発することに決めました。
EzXML.jlの売り
EzXML.jlの売りは、その名前の通り簡単に(easy)にXMLを扱えるようにすることです。そのため、APIや型の設計をもう一度考え直し、かつlibxml2が提供する便利な機能を使えるようにしています。次のコードは、EzXML.jlを使ってXMLのツリーを辿ったり、XPathを使うデモです。
using EzXML # Parse an XML string # (use `readxml(<filename>)` to read a document from a file). doc = parsexml(""" <primates> <genus name="Homo"> <species name="sapiens">Human</species> </genus> <genus name="Pan"> <species name="paniscus">Bonobo</species> <species name="troglodytes">Chimpanzee</species> </genus> </primates> """) # Get the root element from `doc`. primates = root(doc) # Iterate over child elements. for genus in eachelement(primates) # Get an attribute value by name. genus_name = genus["name"] println("- ", genus_name) for species in eachelement(genus) # Get the content within an element. species_name = content(species) println(" └ ", species["name"], " (", species_name, ")") end end println() # Find texts using XPath query. for species_name in content.(find(primates, "//species/text()")) println("- ", species_name) end
EzXML.jlがXMLファイルを読み込むときの関数名はreadxml
です。これは、標準ライブラリにもあるreadcsv
に沿って命名したものです。read
という関数名はJuliaの標準ライブラリではファイルから読み込む際に使われる一般的な動詞ですので、その点でも一貫性があります。文字列からデータを読み込む場合にはparsexml
という関数名になっていて、これもJuliaではparse
が文字列に対して頻繁に使われることを反映されています。ある要素の属性を取るには辞書(Dict
)と同じようにelement[name]
という構文が使えます。値をセットするにはもちろんelement[name] = value
と書きます。他の関数に関しても、一貫性があって一度理解すればあまりレファレンスを引かなくてもプログラミングできるようになっています。
型の設計に関してもなるべくシンプルになるようにしています。EzXML.jlが提供している型は主にEzXML.Document
とEzXML.Node
の2種類だけです。EzXML.Document
はXML文書全体に対応する型で、EzXML.Node
はその構成要素であるノードを表す型です。EzXML.Node
はXML中のノードっぽいものすべてをラップしますので、要素も、テキストも、コメントも、属性もすべてがEzXML.Node
です。そのノードがどのような種類のノードか知るにはnodetype
関数を呼びます。これらの型の名前はコード中で使われることはあまりないですから、パッケージとしてはexportしていません。しかし、必要であれば上のようにモジュール名を付けて使うこともできます。
機能もLightXML.jlより豊富です。名前空間を扱うこともできますし、XPathを使って目的の要素を検索することもできます。XPathはfind
という関数をオーバーロードしていて、find(node, "//foo")
とすればnode
ノード以下のfoo
要素をすべて取得することができます。また、メモリに乗り切らない大規模なXMLファイルを高速にパースするために、streaming readerも提供しています。
最後に、EzXML.jlはメモリリソースを自動的に解放してくれます。普通にJuliaを使っている人はメモリリソースの解放なんて気にしたくないでしょうから、自動的にやってくれれば大変楽です。なので、EzXML.jlではこの辺の面倒を見てあげることにしています。この次では、この実装がどのようになっているのかを説明していきます。
内部実装
まずはJuliaでCのライブラリーをラップする際の基本的なことを確認しておきましょう。JuliaのオブジェクトとC言語の構造体ではメモリレイアウトに互換性があります。例えば、struct point_s { int kind; double x; double y; }
という構造体があれば、Julia側からはimmutable Point; kind::Cint; x::Float64; y::Float64; end
のような型のオブジェクトを定義してやれば、そのメモリレイアウトをそのままJulia側に写すことができます。Cの関数からpoint_s*
のようなポインタ値が返された場合は、Julia側でunsafe_load(ptr::Ptr{Point})
を使ってJuliaのPoint
オブジェクトに変換できます。また、C言語の関数の呼び出しはccall
で実現できます。
EzXML.jlでは、libxml2のノードを扱う時は以下の_Node
型に写しています。これは、libxml2のどんな種類のノードも共通して持っているフィールドです。後で述べますが、最初の_private
フィールドが重要です。
immutable _Node _private::Ptr{Void} typ::Cint name::Cstring children::Ptr{_Node} last::Ptr{_Node} parent::Ptr{_Node} next::Ptr{_Node} prev::Ptr{_Node} doc::Ptr{_Node} end
Juliaはオブジェクトが解放されるときに関数を呼び出す仕組み(finalizer
関数)を提供しており、これを使うことでリソースの解放ができます。しかしこれは、想像するよりちょっと複雑です。何故なら、XMLのノードはlibxml2内で繋がっていますから、誤ってつながっている一部のノードを勝手に解放してしまうと、一貫性が保たれなくなってしまいます。以下のコードを見てみましょう。
function extract_first_child(filename) doc = readxml(filename) r = root(doc) c1 = firstchild(r) return c1 end # foo.xml: # <root> # <child1/> # <child2/> # </root> child = extract_first_child("foo.xml")
この関数ではXML文書中の最初の要素を返しています。関数が返った後はJuliaレベルではdoc
への参照が無くなりますから、Juliaのガーベッジコレクション(GC)はdoc
を解放しようとします。このときにlibxml2のメモリーも解放して良さそうに思えますが、実際にはそうではありません。doc
が参照しているlibxml2の構造体はr
を通してc1
につながっていて、c1
もdoc
へのポインターを持っています。ですので、領域を解放するには、XMLのツリーのどこにもJuliaからの参照が残っていない状態である必要があります。
この問題を回避するために、EzXML.jlではメモリーの解放は常に特定のノードから起きるようにしています。実装上は、Node
型のオブジェクトはptr
というlibxml2のノードを表す構造体へのポインターと自身を管理しているowner
という2つのフィールドを持っています。このowner
オブジェクトが他のノードのメモリー管理を担っているオブジェクトです。
type Node ptr::Ptr{_Node} owner::Node end
Node
のコンストラクターの内部では、以下の様なコードを使って最上位のNode
のポインターを手繰り寄せて、そのノードのowner
に指定しています。つまり、オーナーはXMLツリーの最上位にあるノードがその役割を果たします。
# determine the owner of this node owner_ptr = ptr while unsafe_load(owner_ptr).parent != C_NULL owner_ptr = unsafe_load(owner_ptr).parent end
このowner
フィールドを持つ理由は他にもあります。owner
を持つことで、GCがその子(子孫)ノードより先にオーナーを解放しようとするのを防ぐことができます。上の関数の例で言えば、c1
がdoc
への参照をJuliaのレベルで保持しているため、c1
がある限りdoc
が先に解放されることはありません。さらに途中に挟まれているr
は自身や他のノードのオーナーではないので、Juliaはr
をGCするときにlibxml2で確保されたメモリー領域を解放しません。こうして、常にノードは最上位のノードへの参照を持ち、最上位のノードがメモリーを解放することにして、誤ってGCされるのを防いでいます。
しかし実は、もうひとつ考えなければいけない問題があります。それは、Julia側からあるlibxml2のノードに対して複数の参照を作ってしまうことです。以下の様な簡単なコードを考えてみると、すぐに問題が分かります。
doc1 = readxml("foo.xml")
doc2 = document(root(doc1))
このとき、単純に実装するとdoc1
とdoc2
はlibxml2の同じノードを指しているのに別々のJuliaオブジェクトになっています。すると、あるノードの対して2つのオーナーが存在することになり、両者が最終的に同じノードを解放しようとして後のほうが不正な操作になります。
この問題はlibxml2のノードの構造体にJuliaオブジェクトへのポインターを保持することで解決しています。既に述べたように、libxml2のノードの構造体には_private
と呼ばれるフィールドがあり、ユーザーが自由に使うことができます。そこで、このフィールドにJuliaのNode
オブジェクトに対するポインターを保持しておき、あるノードに対するJulia側のNode
オブジェクトを作る際には、もしあればここからオブジェクトを取り出しています。つまり、既にJulia側であるノードに対するオブジェクトを作成していればそれを使い、無ければ新しく作ってポインターをノードの保持しておくという操作をNode
のコンストラクター内で行っています。こうすることで、あるlibxml2のノードに対するJulia側のNode
オブジェクトは常に高々1個になり、二重に領域を解放してしまうことを防ぐことができます。
最終的に、Node
のコンストラクターは以下の様な感じになります(実際のコードを少し簡略化しています)。
CのポインターからJuliaのオブジェクトを取り出すのがunsafe_pointer_to_objref
で、Juliaのオブジェクトのポインターを取得するのがpointer_from_objref
です。
function Node(ptr::Ptr{_Node}) # return a preallocated proxy object if any str = unsafe_load(ptr) if str._private != C_NULL # found a valid proxy return unsafe_pointer_to_objref(str._private)::Node end # find the owner of this node owner_ptr = ptr while unsafe_load(owner_ptr).parent != C_NULL owner_ptr = unsafe_load(owner_ptr).parent end if ptr == owner_ptr # manage itself node = new(ptr) node.owner = node else # delegate management to its owner owner = Node(owner_ptr) node = new(ptr, owner) end # register finalizer and store the pointer to a node object finalizer(node, finalize_node) unsafe_store!(convert(Ptr{UInt}, ptr), convert(UInt, pointer_from_objref(node))) return node end
finalize_node
は自身がオーナーである場合は、それ以下の管理するノードにつながっているJuliaのオブジェクトを切り離し、libxml2の関数を呼んでノードを解放しています。もしオーナーでなければ、_private
にNULL
ポインターを代入して、JuliaのオブジェクトがGCされてもう存在しないことを示しています。
# Finalize a Node object. function finalize_node(node) node_ptr = node.ptr if node === node.owner # detach pointers to C structs of descendant nodes traverse_tree(node_ptr) do ptr str = unsafe_load(ptr) if has_proxy(str) # detach! unsafe_extract_proxy(str).ptr = C_NULL end end # free the descendants if unsafe_load(node_ptr).typ == DOCUMENT_NODE ccall((:xmlFreeDoc, libxml2), Void, (Ptr{Void},), node_ptr) else ccall((:xmlFreeNode, libxml2), Void, (Ptr{Void},), node_ptr) end elseif node_ptr != C_NULL # indicate the proxy does not exit anymore store_proxy_pointer!(node, C_NULL) end return nothing end
以上で、EzXML.jlの大まかな内部実装がお分かりになるかと思います。実際には、サブツリーを別のツリーへと動かす際にはオーナーをアップデートする必要があるなど細かい処理が多少あるのですが、そのへんまで興味ある方は元のコードを読んでみて下さい。まだ新しいパッケージですので、使ってみてバグや提案などがある場合は報告して下さい。
RユーザーのためのJulia100問100答
R Adevnt Calendar 8日目の記事です。大幅に遅れて大変申し訳ないです。
この記事ではR言語ユーザーのために100問100答形式でJuliaを紹介していこうと思います。
Julia言語
Juliaってどういう言語なの?
Juliaは高レベルでハイパフォーマンスな技術計算のための動的言語だよ。書きやすさと実行速度の両立がウリの言語だよ。
誰が作ってるの?
主にボストンのMITの人達が中心に作っている言語だよ。特にJeff Bezonson, Stefan Karpinski, Viral Shah, Alan Edelmanの4人が初期の重要人物だよ。
自由に使えるの?
Juliaの処理系はMITライセンスで配布されているから、商用でもなんでもかなり自由に使えるよ。
どれくらい速いの?
すごく速いよ!大体C言語の2倍以内くらいの収まる速度だよ。
Rと比較してどうなの?
数倍から数百倍くらい高速かな。特にループや関数呼び出しがたくさんある場合には顕著な差が出るよ。
公式サイトにあるベンチマークの抜粋を載せておくよ (C言語=1.0)。
実例はある?
forループがあるような場合はかなりはっきり差が出るよ。
例えば1からNまでのコラッツ予想を実証する関数collatz
をJuliaとRで書いてみたよ。
function collatz(N) for n in 1:N while n != 1 if iseven(n) n = div(n, 2) else n = 3n + 1 end end end end
collatz <- function (N) { for (n in 1:N) { while (n != 1) { if (n %% 2 == 0) { n = n %/% 2 } else { n = 3 * n + 1 } } } }
N=100,000で試してみるとJuliaはRの300倍ぐらい高速だったよ。
julia> @time collatz(100000)
0.026903 seconds (4 allocations: 160 bytes)
R> system.time(collatz(100000))
user system elapsed
8.429 0.028 8.479
どうして速いの?
LLVMというコンパイラの基盤ライブラリをつかって実行時にコンパイルを行っているからだよ。あとJulia言語自体が最初からパフォーマンスを考慮して設計されているよ。詳しくはBezansonら論文を参照してね。
どこが技術計算に向いてるの?
実行速度が速く、かつ動的プログラミング言語で気軽に実行できるから短いスクリプトで計算するのに向いてるよ。あと、数値計算でよく用いられるデータ型や関数が最初から用意されているから、インストールしてすぐ使い始められるよ。例えば、線形代数のためにBLASやLAPACKの関数が標準で入ってるよ。
Juliaのレポジトリはどこ?
GitHubのレポジトリがここにあるよ: https://github.com/JuliaLang/julia
インストール用のバイナリは?
公式サイトのダウンロードページにWindows, Mac, Linux用のが用意してあるよ: http://julialang.org/downloads/
Juliaのマニュアルやチュートリアルは?
Juliaはドキュメントがよく書かれているから、公式マニュアルを読むのがオススメだよ。 日本語だと私が書いたJuliaのチュートリアルやM.Hiroiさんのサイトがあるよ。
環境
どうやって実行するの?
次のスクリプトをhello.jlとしてファイルに保存したとしよう。
println("hello, world")
Juliaをインストールしてあるならjulia
コマンドがあるはずなので、これにさっきのファイルを渡せばRscript
みたいに実行できるよ。
$ julia hello.jl
hello, world
対話セッション(REPL)を使うには、引数無しでjulia
を実行しよう。そうすると、下のようにRみたいなJuliaのREPLが立ち上がるよ!
$ julia
_
_ _ _(_)_ | A fresh approach to technical computing
(_) | (_) (_) | Documentation: http://docs.julialang.org
_ _ _| |_ __ _ | Type "?help" for help.
| | | | | | |/ _` | |
| | |_| | | | (_| | | Version 0.5.0 (2016-09-19 18:14 UTC)
_/ |\__'_|_|_|\__'_| |
|__/ | x86_64-apple-darwin14.5.0
julia> 1 + 2
3
julia>
ヘルプの調べ方は?
RみたいにREPLに?
と打ってみよう。そうするとプロンプトがhelp?>
に変わるから、知りたい関数の名前をそこに打ってEnterを押そう!
help?> sum
search: sum sum! sumabs summary sumabs2 sumabs! sum_kbn sumabs2! cumsum cumsum! consume cumsum_kbn isnumber isalnum SlotNumber
sum(itr)
Returns the sum of all elements in a collection.
sum(A, dims)
Sum elements of an array over the given dimensions.
sum(f, itr)
Sum the results of calling function f on each element of itr.
キーワード検索をする場合は、help?>
プロンプトでキーワードを"..."
で囲んで入力しよう!
help?> "variance"
Base.var
Base.cov
Base.varm
RStudioみたいなのはあるの?
RStudioほどじゃないかもしれないけど、JunoというAtomエディターをベースにした開発環境が人気だよ。他にはEclipse向けのプラグイン(https://github.com/JuliaComputing/JuliaDT)やVS Code向けのプラグイン(https://github.com/JuliaEditorSupport/julia-vscode)もあるよ。
EmacsやVimの人は?
Emacs用プラグイン(https://github.com/JuliaEditorSupport/julia-emacs)とVim用プラグイン(https://github.com/JuliaEditorSupport/julia-vim)ももちろんあるよ。
ノートブックはあるの?
IJulia.jlを使って、Jupyter Notebookを利用できるよ。すぐ試してみるのにJuliaBoxを使うといいんじゃないかな?
CRANみたいなパッケージのレポジトリはあるの?
あるよ!パッケージの一覧がhttp://pkg.julialang.org/から閲覧できるよ。
パッケージのインストールはどうするの?
Rのinstall.packages
に当たるのがPkg.add
だよ。もちろん依存解決もしてくれるよ。
GitHubにしかないパッケージのインストールは?
Pkg.clone
にレポジトリのURLを入れると野良パッケージでもインストールできるよ。ここでも依存解決をしてくれるよ。
Rを呼び出したいんだけど?
RCall.jlを使おう! RCall.jlはJuliaのコードにRを埋め込んだりRのREPLを使ったりできるよ。
Pythonも呼び出したいんだけど?
PyCall.jlを使おう! ついでに、C++が呼び出せるCxx.jlとかMATLABが呼び出せるMATLAB.jlなんてのもあるよ。
CやFortranのコードの呼び出しは?
JuliaはCやFortranのコードを呼び出す仕組みを標準で用意しているよ。 詳しくは後で説明するよ。
Rのサンプルデータを読みたいんだけど?
RDatasets.jlにRでよく使われるデータセットがパッケージされているよ!
ワークスペースを保存するにはどうするの?
JLD.jlの@save
マクロを使おう。
julia> using JLD
julia> n = 5
5
julia> x = randn(n, n)
5×5 Array{Float64,2}:
-2.18618 0.0831716 2.19386 0.806268 1.11444
-0.99689 -0.187922 -0.138358 -0.141601 0.19058
0.0971969 0.149858 -1.19826 1.15507 1.0969
-0.140956 -0.727159 0.780267 1.17143 0.979918
-0.161376 -0.162103 -0.175158 -0.402853 0.916248
julia> @save "workspace.jld"
ワークスペースを復元するには?
同じくJLD.jlの@load
マクロを使おう。
julia> using JLD
julia> @load "workspace.jld"
3-element Array{Symbol,1}:
:ans
:n
:x
julia> n
5
julia> x
5×5 Array{Float64,2}:
-2.18618 0.0831716 2.19386 0.806268 1.11444
-0.99689 -0.187922 -0.138358 -0.141601 0.19058
0.0971969 0.149858 -1.19826 1.15507 1.0969
-0.140956 -0.727159 0.780267 1.17143 0.979918
-0.161376 -0.162103 -0.175158 -0.402853 0.916248
コミュニティ
Juliaのコミュニティはどんな感じなの?
Juliaのコミュニティは小さいけれど、とてもオープンなコミュニティを形成しているよ。 特に初心者に優しい(Juliaユーザはだいたい初心者)ので、気軽にコミュニケーションを取れるよ。
質問はどこですればいいの?
英語ならDiscourseかStack Overflowだね。Stack Overflowで質問する時は"julia-lang"のタグをつけるようにしよう。 Discourseの方がJuliaハッカーが見てる率が高い気がするから、Juliaに限定した質問ならDiscourseの方がオススメかも。
日本語話者のコミュニティは?
日本ではJulia Tokyoが知る限り唯一のコミュニティかな? Facebook・Slack・GitHubのアカウントがあるし、不定期でミートアップも開催してるよ。
他にはどんなコミュニティがあるの?
Juliaのコミュニティは各ドメインに分かれてゆるくつながっているよ。 コミュニティのリストはここ(http://julialang.org/community/)にあるよ。 いくつか挙げると以下のドメインのコミュニティは結構活発かな。
ドメイン | GitHub |
---|---|
統計 | JuliaStats |
最適化 | JuliaOpt |
データベース | JuliaDB |
GPU | JuliaGPU |
生物学 | BioJulia |
微分 | JuliaDiff |
グラフィクス | JuliaGraphics |
計量経済学 | QuantEcon |
これらで活発に開発されているパッケージは概ね安心して使えるかな。
どんな人達がJuliaを使ってるの?
Juliaの使い方は色々あるけれど、やっぱり数値計算目的に使っているユーザが多いかな。 Julia Computing, Inc.の事例を見ると、金融・経済分野を始めとしたデータ分析界隈で使われ始めている印象があるな。
HPCでの事例や科学論文でJuliaを使った計算の実装も最近よく見るようになってきたよ。
他には、StanfordやMITでは技術計算の授業にもJuliaが使われているよ。
トークとかの映像は?
YouTubeのチャンネルにカンファレンスの動画などがたくさんあるよ!
https://www.youtube.com/user/JuliaLanguage
Julia界のHadley Wickhamは?
Julia界には(良くも悪くも)あれくらい目立っている人はいないかな。 Juliaの重要なパッケージは各organizationがコミュニティとして開発・管理してるよ。
文法
Rの構文との対応関係を教えて?
いいよ!
代入
# R x <- 100
# Julia x = 100
分岐
# R if (x == 0) { print("zero") } else if (x < 0) { print("negative") } else { print("positive") }
# Julia if x == 0 println("zero") elseif x < 0 println("negative") else println("positive") end
for文
# R for (i in 1:10) { print(i) if (i > 5) { break } }
# Julia for i in 1:10 println(i) if i > 5 break end end
関数
# R add <- function (x, y) x + y add <- function (x, y) { return (x + y) }
# Julia add(x, y) = x + y function add(x, y) return x + y end
ライブラリ読み込み
# R library(ggplot2)
# Julia using DataFrames
たまにJuliaのコードにある@
は何?
@
から始まるコードはJuliaのマクロ呼出しだよ。
マクロはJuliaのコードを別のコードに書き換えたりするメタプログラミングの一種だよ。
例えば次のコードに出てくる@inbounds
は配列アクセスの境界チェックを無くして少し高速化することができるし、@show
はその時の値をいい感じに表示してくれるよ。
x = randn(10) s = 0 @inbounds for i in 1:endof(x) s += x[i] end @show s
たまにJuliaのコードにあるr"..."
は何?
それもJuliaのマクロだよ!
標準ライブラリではr"..."
と書くことで文字列でなく正規表現を作ることができるよ。
julia> r"foo"
r"foo"
julia> typeof(r"foo")
Regex
julia> b"foo"
3-element Array{UInt8,1}:
0x66
0x6f
0x6f
julia> typeof(b"foo")
Array{UInt8,1}
この仕組はユーザからも拡張可能なので、例えばBio.jlではDNA配列を作ったりするのに使ってるよ。
julia> using Bio.Seq
julia> dna"ACGTAG"
6nt DNA Sequence:
ACGTAG
データ
主な数値型は?
整数はInt
型、倍精度浮動小数点数はFloat64
型がデフォルトでよく使われるよ。
Int
は32-bit環境なら32-bit、64-bit環境なら64-bitで表現されるよ。
Juliaは数値の型も豊富で、8-bitから128-bitまで符号の有り無しの組み合わせがすべて用意されているし、複素数(Complex{T}
)や有理数(Rational{T}
)もあるよ。
TRUE
やFALSE
は?
Juliaでは全部小文字のtrue
とfalse
だよ。
欠損値は扱えるの?
JuliaではNullable{T}
という型があって、型T
の欠損値を表すよ。
nothing
という値もあるけどパフォーマンス上の理由であまりお勧めしないよ。
JuliaもRみたいに常に配列として数値を扱うの?
Juliaでは単一の値と配列は区別されるよ。例えば、3
と書いたら整数の3を意味するし、[3]
と書いたら1つの整数3からなるベクトルだよ。
3
と書いたら浮動小数点数じゃなくて整数なの?
そうだよ、Rと違ってJuliaでは整数として扱われるよ。
じゃぁ浮動小数点数の3を作るには?
3.0
と書けば浮動小数点数として扱われるよ。
どうやって確認するの?
Rみたいにtypeof
という関数を使おう!
julia> typeof(3)
Int64
julia> typeof(3.0)
Float64
R> typeof(3)
[1] "double"
R> typeof(3.0)
[1] "double"
演算子はどんな感じ?
ほとんどRと同じだよ。
julia> 3 + 4
7
julia> 3 - 4
-1
julia> 3 * 4
12
julia> 3 / 4
0.75
R> 3 + 4
[1] 7
R> 3 - 4
[1] -1
R> 3 * 4
[1] 12
R> 3 / 4
[1] 0.75
Bool演算(&&
, ||
, !
)も同じだよ。排他的論理和はv0.5では$
だけど次期バージョンでxor
という名前に変更される予定だよ。
list
はあるの?
Juliaで一番近いのはDict
かな?
Pythonのdict
と同じハッシュを使った連想配列だよ。
julia> d = Dict("one" => 1, "two" => 2, "three" => 3)
Dict{String,Int64} with 3 entries:
"two" => 2
"one" => 1
"three" => 3
julia> d["two"]
2
data.frame
は?
残念ながらJuliaの標準ライブラリにはないよ。でもDataFrames.jlというパッケージがJuliaのデータフレームの標準的な実装だよ。 この辺は後で詳しく述べるよ。
多次元行列も扱えるの?
もちろん扱えるよ! Rのmatrix
みたいなことが次のようにできるよ。
julia> [1 2 3
4 5 6]
2×3 Array{Int64,2}:
1 2 3
4 5 6
julia> [1 2 3; 4 5 6] # セミコロンは改行と同じ扱い
2×3 Array{Int64,2}:
1 2 3
4 5 6
R> matrix(c(1, 2, 3, 4, 5, 6), nrow=2)
[,1] [,2] [,3]
[1,] 1 3 5
[2,] 2 4 6
空の配列の作成は?
ベクトル(1次元配列)なら[]
だよ。
julia> []
0-element Array{Any,1}
ただ要素の型が分かっている場合はそれを教えてあげるとJuliaは高速に動くよ。
julia> Int[]
0-element Array{Int64,1}
julia> Float64[]
0-element Array{Float64,1}
matrix(0, 3, 4)
みたいに0で初期化された配列の作成は?
zeros
関数を使おう!
julia> zeros(3, 4) # 3行4列の行列
3×4 Array{Float64,2}:
0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0
julia> zeros(Int, 3, 4)
3×4 Array{Int64,2}:
0 0 0 0
0 0 0 0
0 0 0 0
1で初期化された配列ならones
だよ。
配列要素へのアクセスはどうするの?
基本的にRと同じだよ。Juliaの配列も1始まりでcolumn majorだよ。
ただ、Rだとx[i,]
とかx[,j]
と書くところをJuliaではx[i,:]
とかx[:,j]
とか書くよ。
julia> x = [1 2 3; 4 5 6]
2×3 Array{Int64,2}:
1 2 3
4 5 6
julia> x[1,2]
2
julia> x[2,2:3]
2-element Array{Int64,1}:
5
6
julia> x[:,2]
2-element Array{Int64,1}:
2
5
2:3
みたいなのはベクトル?
Juliaではm:n
はm
からn
までの範囲(Range
)を表すよ。
julia> 2:3
2:3
julia> typeof(2:3)
UnitRange{Int64}
ベクトルと違ってメモリーを消費しないので、長い範囲でも一瞬で作れるよ。
character
みたいな文字列型は?
JuliaではString
型があるよ。
文字列はUnicodeも使える?
使えるよ! JuliaではUTF-8でエンコーディングしてあるよ。
文字列操作関数を教えて?
stringrパッケージの関数と大雑把に対応を付けたよ。
R (stringr) | Julia |
---|---|
str_c |
string |
str_detect |
ismatch |
str_extract |
match |
str_dup |
^ |
str_split |
split |
str_replace |
replace |
str_trim |
strip , lstrip , rstrip |
str_pad |
lpad , rpad |
str_to_lower |
lowercase |
str_to_upper |
uppercase |
as.double
とかas.integer
みたいな型変換はどうするの?
Juliaの型変換は基本的にconvert
関数を使うよ。
julia> convert(Float64, 42)
42.0
julia> convert(Int64, 42.0)
42
ただし文字列をパースして数値にするにはparse
を使うよ。
julia> parse(Int, "42")
42
julia> parse(Float64, "42")
42.0
行名や列名はどうつけるの?
Juliaの標準では行名や列名をサポートしていないよ。 でもNamedArrays.jlという配列に名前を付けられるようにするパッケージはあるよ。
欠損値のある配列はどう扱えばいいの?
欠損値の扱いはまだJulia界隈でも完全な合意に達してない課題なので、ちょっとむずかしいな。
前にも述べたようにNullable{T}
というデータ型があるけど、これに関する操作が将来変わる可能性があるよ。
でも、基本的にはNullableArrays.jlを使おう!
因子(factor)を扱うにはどうするの?
これもまだJuliaではちょっと扱いが難しいけど、CategoricalArrays.jlというパッケージがあるよ。 NullableArrays.jlやCategoricalArrays.jlは近いうちにDataFrames.jlに採用される予定のようだよ。
関数
関数はどうやって定義するの?
function
キーワードを使ったスタイルと短縮記法があるよ。
# 通常の記法 function add(x, y) return x + y end # 短縮記法 add(x, y) = x + y
デフォルト引数やキーワード引数は?
あるよ! でもJuliaではデフォルト引数とキーワード引数で記法が分かれているよ。
# デフォルト引数 function default(x, y=1, z=10) return x + y + z end # キーワード引数 (セミコロンに注目) function keyword(x; y=1, z=10) return x + y + z end
試しに実行してみよう
julia> default(1)
12
julia> default(1, 0)
11
julia> default(1, 0, 5)
6
julia> keyword(1)
12
julia> keyword(1, z=5)
7
julia> keyword(1, z=5, y=0)
6
クロージャも作れる?
もちろん!
function counter() i = 0 function x(n=1) i += n return i end end
julia> c = counter()
(::x) (generic function with 2 methods)
julia> c()
1
julia> c()
2
julia> c()
3
julia> c(3)
6
julia> c()
7
+
とかは関数じゃないの?
関数だよ! 前置記法で書くこともできるよ!
julia> +(1, 2)
3
R> `+`(1, 2)
[1] 3
じゃぁRみたいにif
とかfor
も関数なの?
Juliaではif
・for
・while
は関数じゃないよ。
遅延評価はあるの?
Juliaには遅延評価はなくてすべて正格評価だよ。つまり関数に渡された引数は呼び出し前にかならず評価されるよ。
apply
みたいな関数を配列の要素に適用する関数は?
map
を使おう! ただし関数が最初の引数だよ。
julia> map(log, [1,2,3])
3-element Array{Float64,1}:
0.0
0.693147
1.09861
julia> map(exp, [1,2,3])
3-element Array{Float64,1}:
2.71828
7.38906
20.0855
R> sapply(c(1,2,3), log)
[1] 0.0000000 0.6931472 1.0986123
R> sapply(c(1,2,3), exp)
[1] 2.718282 7.389056 20.085537
ただ最近はブロードキャストを使う方法が一般的なので、こちらを使ったほうが良いよ。
julia> log.([1,2,3])
3-element Array{Float64,1}:
0.0
0.693147
1.09861
julia> exp.([1,2,3])
3-element Array{Float64,1}:
2.71828
7.38906
20.0855
この関数名の後のドットは何?
ブロードキャスト関数呼び出しの構文糖衣だよ。例えばf.(x)
はbroadcast(x, f)
に変換されるよ。
broadcast
はmap
の凄い版だよ。
ベクトル同士の足し引きとかは?
+
や-
がそのまま使えるよ。でも.
を付けた.+
や.-
なら左右で配列のサイズが違った場合でもRのように動いてくれるよ。
julia> [1,2,3] + [4,5,6]
3-element Array{Int64,1}:
5
7
9
julia> [1,2,3] .+ [4]
3-element Array{Int64,1}:
5
6
7
julia> [1,2,3] .+ 4
3-element Array{Int64,1}:
5
6
7
R> c(1, 2, 3) + c(4, 5, 6)
[1] 5 7 9
R> c(1, 2, 3) + c(4)
[1] 5 6 7
R> c(1, 2, 3) + 4
[1] 5 6 7
書いた関数が思ったより遅いんだけど?
一番よくある原因が関数内で使われている変数の型が不安定になっているせいだよ。
例えば次のコードはxs
の要素がInt
のときには高速に動くけど、浮動小数点数の配列とかを渡すとすごく遅いよ。
function sumup(xs) s = 0 for x in xs s += x end return s end
julia> xs = rand(1:10, 10000);
julia> sumup(xs); # ウォームアップ
julia> @time sumup(xs)
0.000011 seconds (5 allocations: 176 bytes)
54649
julia> xs = rand(10000);
julia> sumup(xs); # ウォームアップ
julia> @time sumup(xs)
0.000337 seconds (30.00 k allocations: 468.906 KB)
5015.030044010504
原因は変数s
は整数型で初期化されてるのに、ループの中で浮動小数点数と足すことでs
の型が整数か浮動小数点数かJuliaには分からなくなってしまうせいだよ。これをJulia界隈では型の不安定さ(type instability)と言ってるよ。
これを回避するにはeltype
とzero
関数を次のように使うことで、s
の型をxs
に合わせて設定できるよ。
function sumup(xs) s = zero(eltype(xs)) for x in xs s += x end return s end
こうすることでさっきより何十倍も速くなったね!
julia> sumup(xs);
julia> @time sumup(xs)
0.000015 seconds (5 allocations: 176 bytes)
5015.030044010504
関数のプロファイルを取るにはどうするの?
@profile
マクロと、Profile.print()
をあわせて使おう!
@profile <code>
でコード実行のプロファイルをとって、Profile.print()
でプロファイル結果をプリントできるよ。
C言語の関数の呼び出し方法を教えて!
Juliaにはccall
というC言語を呼び出す専用の命令があるよ。例えば、GNU libcに実装されているexp
関数を呼び出すなら次のように書けるよ。
julia> ccall((:exp, "libc"), Float64, (Float64,), 1.0)
2.718281828459045
julia> ccall((:exp, "libc"), Float64, (Float64,), -1.0)
0.36787944117144233
ccall
は最初に呼び出す関数を指定して、その後に返り値の型と引数の型を指定して、最後に引数を渡すよ。上の例で言えば、(:exp, "libc")
が呼び出す関数で、Float64
が返り値の型、(Float64,)
が引数の型で、1.0
や-1.0
がCの関数に渡される実際の引数だよ。
Juliaは数値以外にも文字列やより複雑な構造体をCのライブラリとやり取りすることができるように設計されているので、より詳しくはマニュアルのここを読んでみてね。 具体的な例としてはJuliaの標準ライブラリやLibz.jlなんかが勉強になるよ。
型システムとメソッド
Juliaの型システムはどうなってるの?
Juliaの型は大きく具体型(concrete type)と抽象型(abstract type)に分けられていて、メソッドの選択は多重ディスパッチというシステムを中心に構築されているよ。
具体型とか抽象型って何?
具体型はインスタンス化(オブジェクトを作れる)できる型で、抽象型はそうでない型だよ。
例えばFloat64
は1.0
など値を作れるけどReal
という抽象型は値を作れないよ。
値を作れない抽象型は何の役に立つの?
複数の型をまとめるのに使われるよ。KeitaNakamuraさんのコードを使わせてもらうと、例えばReal
を頂点とする型は次のような構造になっているよ。
julia> print_tree(Real)
Real
├─── AbstractFloat
| ├─── BigFloat
| ├─── Float16
| ├─── Float32
| └─── Float64
├─── Integer
| ├─── BigInt
| ├─── Bool
| ├─── Signed
| | ├─── Int128
| | ├─── Int16
| | ├─── Int32
| | ├─── Int64
| | └─── Int8
| └─── Unsigned
| ├─── UInt128
| ├─── UInt16
| ├─── UInt32
| ├─── UInt64
| └─── UInt8
├─── Irrational{sym}
└─── Rational{T<:Integer}
function print_tree(typ) println(typ) print_tree(typ, []) end function print_tree(typ, level) stypes = subtypes(typ) for stype in stypes if stype !== stypes[end] println(join(level) * "├─── $stype") push!(level, "| ") else println(join(level) * "└─── $stype") push!(level, " ") end print_tree(stype, level) pop!(level) end end
クラス(型)の定義はどうやるの?
RのS4に似た仕組みがJuliaにはあるよ。例えば以下の2つの定義は大体対応しているよ。
# R setClass("Person", representation(name = "character", age = "numeric"))
# Julia type Person name::String age::Int end
Juliaはデフォルトコンストラクタを準備しているから、オブジェクトは次のように作成できるよ。
julia> hadley = Person("Hadley", 31)
Person("Hadley",31)
julia> hadley.name
"Hadley"
julia> hadley.age
31
ただし、Juliaのtype
で宣言したオブジェクトは変更可能なのでその振る舞いはRの参照クラス(RC)に近いかな。
継承は?
JuliaはsetClass
のcontains
引数に当たるような継承の仕組みはないよ。つまり、複数の型間で構造を継承する方法はないよ。
その代わり、定義する型をある抽象型のサブタイプにする仕組みはあって、Juliaではこちらを使うよ。
例えば、Juliaで要素の数が3つのベクトルっぽい型を作る場合には、次のようにAbstractVector
のサブタイプとして宣言すると良いよ。
type Triplet{T} <: AbstractVector{T} x::T y::T z::T end
ここで、T
と書かれているのは型パラメータで、T = Int64
やT = Float64
が代入されて具体的な型になるよ。
また、各抽象型には実装していることが期待されるメソッドがあるので、それを実装するべきだよ (AbstractVector
ならBase.size
・Base.getindex
などを実装すべき)。
メソッドの定義は?
通常の関数定義に型指定したものがRのsetMethod
に近いよ。
# R setGeneric("greet", function(x) standardGeneric("greet")) setMethod("greet", signature("Person"), function(x) paste("Hello, ", x@name))
# Julia function greet(x::Person) return string("Hello, ", x.name) end
型指定は複数に引数に対しても指定できるし、抽象型も指定できるよ。これをJuliaでは多重ディスパッチと言うよ。
前処理
CSV/TSVファイルの読み込みは?
3通りの方法があるよ!
- 標準ライブラリの
readcsv
- DataFrames.jlの
readtable
- CSV.jlの
CSV.read
どれがいいの?
標準ライブラリのreadcsv
は標準なのですぐ使えるけど、行列としてデータを読み込むから、場合によってはちょっと使いづらいかな?
DataFrames.jlのreadtable
はデータフレームを読み込むけど、将来的に使われなくなる可能性が高いよ。
なので、現在ではCSV.jlのCSV.read
を使うのがオススメだよ。
julia> using DataFrames
julia> using CSV
julia> head(CSV.read("iris.csv"))
6×5 DataFrames.DataFrame
│ Row │ SepalLength │ SepalWidth │ PetalLength │ PetalWidth │ Species │
├─────┼─────────────┼────────────┼─────────────┼────────────┼──────────┤
│ 1 │ 5.1 │ 3.5 │ 1.4 │ 0.2 │ "setosa" │
│ 2 │ 4.9 │ 3.0 │ 1.4 │ 0.2 │ "setosa" │
│ 3 │ 4.7 │ 3.2 │ 1.3 │ 0.2 │ "setosa" │
│ 4 │ 4.6 │ 3.1 │ 1.5 │ 0.2 │ "setosa" │
│ 5 │ 5.0 │ 3.6 │ 1.4 │ 0.2 │ "setosa" │
│ 6 │ 5.4 │ 3.9 │ 1.7 │ 0.4 │ "setosa" │
詳しい使い方はCSV.jlのヘルプを参照してね。
SQLiteのデータを読みたいんだけど?
SQLite.jlを使おう。SQLを発行した結果はDataFrame
として返してくれるよ。
JSONを扱うには?
JSON.jlを使おう。
julia> using JSON
julia> JSON.parse("""{"foo": 100, "bar": [1.1, 2.0]}""")
Dict{String,Any} with 2 entries:
"bar" => Any[1.1,2.0]
"foo" => 100
XMLやHTMLのデータを処理するには?
拙作のEzXML.jlがオススメだよ。長く使われているLightXML.jlというのもあるけれど、機能が少ないし、あまりインターフェースがきれいじゃないね...
julia> using EzXML
julia> readxml("ex1.xml")
EzXML.Document(EzXML.Node(<DOCUMENT_NODE@0x00007fcda9030400>))
julia> readxml("ex1.xml") |> print
<?xml version="1.0" encoding="UTF-8"?>
<bookstore>
<book category="COOKING" tag="first">
<title lang="en">Everyday Italian</title>
<author>Giada De Laurentiis</author>
<year>2005</year>
<price>30.00</price>
</book>
<book category="CHILDREN">
<title lang="en">Harry Potter</title>
<author>J K. Rowling</author>
<year>2005</year>
<price>29.99</price>
</book>
</bookstore>
でもHTMLに関してはGoogleのgumboをラップしたGumbo.jlの方が良いかな?
dplyrは?
残念ながらdplyrほどの完成度のパッケージは今のところJuliaには無いかな。 ただJuliaの人にはR使いも多いので、dplyrの利便性は誰もが認めるところだね。 それで最近は似たことをJuliaでやろうとしている人達が出てきて、いくつかパッケージができているよ。
Juliaのデータフレームの近況と将来についてはここにまとまってるよ。
UTF-8以外のファイルを読み込みたいんだけど?
StringEncodings.jlを使おう! libiconvを使ってるのでSJISみたいなアレなエンコーディングでもちゃんと扱えるよ!
統計/機械学習
平均とか分散の計算は?
mean
, median
, var
, cov
, cor
はRと同名の関数が標準ライブラリにあるよ。標準偏差はsd
でなくstd
になってるので注意だよ。
rowSums
とかcolSums
とかは?
sum(array, dim)
が使えるよ。
julia> x
2×3 Array{Int64,2}:
1 2 3
4 5 6
julia> sum(x, 1) # colSums
1×3 Array{Int64,2}:
5 7 9
julia> sum(x, 2) # rowSums
2×1 Array{Int64,2}:
6
15
他にもmean
やvar
なんかでも同じことができるよ。
乱数生成はどうするの?
一様乱数runif(n)
に対応するのはrand(n)
、正規分布rnorm(n)
に対応するのはrandn(n)
だよ。
ガンマ分布とかポアソン分布みたいな他の分布は?
Distributions.jlを使おう。
PCA(主成分分析)をしたいんだけど?
MultivariateStats.jlを使おう。
julia> using MultivariateStats
julia> using RDatasets
julia> iris = dataset("datasets", "iris")
julia> pca = fit(PCA, Matrix(iris[:,1:4])') # 行列に変換して転置
PCA(indim = 4, outdim = 3, principalratio = 0.99479)
julia> transform(pca, Matrix(iris[:,1:4])')
3×150 Array{Float64,2}:
2.68413 2.71414 2.88899 2.74534 2.72872 2.28086 … -1.94411 -1.52717 -1.76435 -1.90094 -1.39019
0.319397 -0.177001 -0.144949 -0.318299 0.326755 0.74133 0.187532 -0.375317 0.0788589 0.116628 -0.282661
0.0279148 0.210464 -0.0179003 -0.0315594 -0.0900792 -0.168678 -0.177825 0.121898 -0.130482 -0.723252 -0.36291
GML(一般化線形モデル)を扱いたいんだけど?
GLM.jlを使おう。Rでも有名なDouglas Bates先生の作ったパッケージだから安心だね!
Lassoは?
Julia実装のLasso.jlかglmnetをラップしたGLMNet.jlが良いと思うよ。
Stanを使いたいんだけど?
Stan.jlを使おう!
行列積の計算は?
Rの%*%
がJuliaでは*
だよ。
julia> A = randn(2, 3);
julia> B = randn(3, 4);
julia> A * B
2×4 Array{Float64,2}:
2.42828 -1.39917 0.28215 -1.61981
0.292669 -0.411146 0.234041 1.82686
R> A <- matrix(rnorm(6), 2)
R> B <- matrix(rnorm(12), 3)
R> A %*% B
[,1] [,2] [,3] [,4]
[1,] 0.8448714 -0.6084743 -0.7458281 -0.653775
[2,] 0.2133118 -0.9305351 0.3933490 1.105863
solve
みたいな線形方程式の解の計算は?
A * X = B
を解くには、A \ B
と計算するよ。
julia> A = randn(3, 3)
3×3 Array{Float64,2}:
0.864382 -0.945314 -0.496754
0.119709 0.83993 -0.962021
-0.737123 0.221293 0.205341
julia> X = randn(3, 2)
3×2 Array{Float64,2}:
-0.224859 1.64374
0.0726233 0.898944
-0.387804 0.697085
julia> B = A * X
3×2 Array{Float64,2}:
-0.0703728 0.224752
0.407157 0.28121
0.102188 -0.869566
julia> A \ B
3×2 Array{Float64,2}:
-0.224859 1.64374
0.0726233 0.898944
-0.387804 0.697085
julia> A \ B ≈ X
true
A
が実対称行列であることが分かってるんだけど?
cholfact(A)
でコレスキー分解してからでも同様に計算できるよ!
julia> A = randn(3, 3);
julia> A = A'A
3×3 Array{Float64,2}:
1.13074 2.65079 -0.0602727
2.65079 7.18991 -0.474305
-0.0602727 -0.474305 0.488933
julia> issymmetric(A)
true
julia> X = randn(3, 2)
3×2 Array{Float64,2}:
0.809623 0.595824
-0.25161 -2.33654
0.97426 1.46868
julia> B = A * X;
julia> A \ B
3×2 Array{Float64,2}:
0.809623 0.595824
-0.25161 -2.33654
0.97426 1.46868
julia> cholfact(A) \ B
3×2 Array{Float64,2}:
0.809623 0.595824
-0.25161 -2.33654
0.97426 1.46868
log(sum(exp(x)))
とかlog(exp(x)+1)
とかを安定的に計算したいんだけど?
StatsFuns.jlに標準ライブラリにはないけど、数値計算でよく見る関数の実装があるからそっちを探してみよう。
julia> x = randn(10)
10-element Array{Float64,1}:
0.0410153
-0.407537
1.26758
-0.52522
-0.0969684
-2.52716
-0.633663
-1.2695
-0.474312
-0.800294
julia> log(sum(exp(x)))
2.165782421890441
julia> logsumexp(x)
2.165782421890441
julia> log(sum(exp(1000x))) # 1000倍すると計算できない
Inf
julia> logsumexp(1000x) # StatsFuns.jlの実装なら大丈夫
1267.5802777924926
プロット作るにはどうするの?
残念ながら、今のところJuliaには決定的なプロットのための仕組みがあるわけではないので、次のものを参考に好みに合ったのを使おう!
- Gadfly.jl
- 最も歴史のあるJuliaの
- 一番スターが多い
- Plots.jl
- 新興で急成長中のライブラリ
- 開発が活発
- PlotlyJS.jl
- Plot.lyのラッパーライブラリ
- インタラクティブな可視化が可能
- PyPlot.jl
- Pythonのmatplotlibのラッパーライブラリ
- 最も高機能かも
- GR.jl
- GRフレームワークのラッパーライブラリ
- パフォーマンスに優れる
- UnicodePlots.jl
- Unicode文字を使ってターミナル上でプロットできるライブラリ
- ごく簡単なプロットをチラッと見たいときに向いてる
関数の最大・最小値を探したいんだけど?
Optim.jlを使おう!
julia> using Optim
julia> rosenbrock(x) = (1.0 - x[1])^2 + 100.0 * (x[2] - x[1]^2)^2
rosenbrock (generic function with 1 method)
julia> optimize(rosenbrock, randn(2))
Results of Optimization Algorithm
* Algorithm: Nelder-Mead
* Starting Point: [0.4289709134439212,-1.3041568238701216]
* Minimizer: [1.0000314591723356,1.0000648302370259]
* Minimum: 1.354834e-09
* Iterations: 69
* Convergence: true
* √(Σ(yᵢ-ȳ)²)/n < NaN: false
* Reached Maximum Number of Iterations: false
* Objective Function Calls: 136
R> rosenbrock <- function (x) (1 - x[1])^2 + 100 * (x[2] - x[1]^2)^2
R> optim(rnorm(2), rosenbrock)
$par
[1] 0.9995882 0.9991591
$value
[1] 2.002849e-07
$counts
function gradient
129 NA
$convergence
[1] 0
$message
NULL
Juliaの機械学習パッケージはどうなってるの?
Juliaユーザには機械学習を常用する人も多いけれど、PythonやRほど充実はしてないよ。 一応、ScikitLearn.jlというscikit-learnのラッパーがJuliaのパッケージになってるので、scikit-learnでできることはできるはずだよ。
Gradient boostingとかは?
XGBoost.jlが使えるよ! 名前から分かる通りXGBoostへのJuliaのインターフェースだよ。
Deep learningのライブラリは?
MXNet.jlというMXNetのパッケージがあるよ。またDMLCだね。他にはMocha.jlというJulia製のフレームワークもあるけれど、メインの開発者がmxnetの方に移ったようなので今はそれほど開発が活発でもないかな。 他にはTensorFlow.jlというPythonのTensorFlowのラッパーライブラリもあるよ。
もっとJuliaっぽいやつだと、Transformations.jlとかNAISTの進藤先生が作っているMerlin.jlというパッケージもあるよ。
結局、Juliaでデータ解析はどうなの?
Juliaは機械学習のアルゴリズムを組むには良い言語だと思うけど、ライブラリやツールもとても少ないね。 なので現状、インタラクティブな探索的データ解析(exploratory data analysis)をするような場合にはJuliaよりRの方をオススメするよ。
ただ、最近JuliaMLという機械学習に特化したコミュニティができてJuliaで気軽に機械学習ができるようにしていたり、データフレームに本格的に見直そうとしていたりとか、改善に向けて動いているよ。今ならやることがたくさんあるので、活動に参加してみてはどうかな?
Julia現状確認 (言語編)
これはJuliaアドベントカレンダー1日目の記事です。2016年におけるJuliaの「現状」を確認していこうと思います。
Julia言語とは
Juliaとは一体どのような言語でしょうか。julialang.orgの最初の説明を用いると、「他の技術計算環境のユーザに馴染みのある構文を備えた、高レベル・高パフォーマンスな技術計算にための動的プログラミング言語」("a high-level, high-performance dynamic programming language for technical computing, with syntax that is familiar to users of other technical computing environments.")です。ここで言う「他の技術計算環境」とは、MATLABやPythonのことを指していると考えられます。実際、Juliaの構文や関数名はMATLABとの類似性が高いため、Octaveのようなオープンソース版のMATLAB代替と捉えられることも多いようですが、プログラミング言語としてはかなり性質が異なります。また、Pythonとは文法や関数名の面ではそれほど類似性は無いように思いますが、似ている点も多数あります。
Juliaの最大の特徴を挙げるならば、動的プログラミング言語でありながらC言語などに匹敵するほど実行が高速であるという点でしょう。多くの技術計算では関数の呼び出しやループが頻繁に用いられますが、Juliaでは実行時コンパイル(JIT)のおかげでこれらの実行コストがネイティブコードにコンパイルするプログラミング言語並になっています。しかし、Juliaが行う実行時のコンパイルは、他の言語でよくみられるJITコンパイラとは性質を異にするようです。PyPy(PythonのJITコンパイラ実装)やLuaJIT(LuaのJITコンパイラ実装)はTracing JITと呼ばれるもので、実行時にプロファイルを行ってホットスポットを同定し、実行時の情報を集め、コンパイルと最適化を行います。しかし、JuliaのJITコンパイラはコードの実行前に型推論と最適化を行い、その後コードを実行します。型推論の結果、型が決定できる部分では強力な最適化を行い、曖昧性が残る部分では実行時の型に合わせて実行できるようにしています。
Julia言語の現状
Julia言語の実装であるJulia処理系はすべてGitHub上で開発されています。レポジトリはこちら( https://github.com/JuliaLang/julia)です。2016年12月1日現在の最新のリリース版はv0.5.0で、9月20日頃リリースされたものです。次期バージョンはv0.6でほぼ間違いないと思います。気になるのはv1.0のリリース時期ですが、今年のJuliaConでJuliaの生みの親のひとりであるStefan Karpinski氏が宣言したことによると、Julia 1.0のリリースは2017年を予定しているようです。次期バージョンのv0.6はv1.0のα版的な意味合いで開発されているようなので、v0.6が出てすぐv1.0がリリースされるものと思われます。しかし実は、v1.0が出てしばらく1.0のままということはなく、1-2年のうちにv2.0も出すつもりのようです。
振り返って、現行のv0.5ではどのような変更があったのかを確認してみます。まず性能面で大きな影響があったのが、クロージャが高速になったことです。v0.4までのJuliaでは、関数がすべて一つの型にまとめられており最適化の恩恵を受けにくかったのですが、v0.5からは関数が別々の型を持つようになりました。これによって、無名関数などの関数の呼び出しが高速化され、安心して使えるようになりました。例えば、map(x -> 2x, array)
のようなコードがforループと同等の速度が出るようになります。また、ASCIIとUTF-8に分かれていた文字列型がString
型に統一され、内部のエンコーディングを気にせず使えるようになりました。これで、Julia内では常にUTF-8のエンコーディングのみを考慮すればよく、プログラムを書くのが大変楽になりました。さらに、新機能として、実験的にマルチスレッディングのサポートが入りました。今のところ、Julia内部でまだマルチスレッディングに対応していない部分(例えばIOなど)があるので、常に問題なく使えるわけではありませんが、競合などがない場合にはforループに@threads
をつけるだけで中身を並列に実行してくれます。これらの変更の他にも、generator式(i for i in 1:10
などが式として受け渡しできる)やfusing構文(sin.(cos.(x))
などがベクトルx
に対して一時配列を作らない)などの新しい機能が導入されています。
では、v0.6はどうなるでしょうか。v0.6はv1.0のリリースのためにいつもより早めにfeature freezeが入るようですので、そろそろ新機能や変更点などが出揃うと思います。その中で、おそらく一番大きな変更が型システムの変更になるでしょう。実装はここ(https://github.com/JuliaLang/julia/pull/18457)にプルリクがあるので、すぐに試すことができます。これを少し詳しく説明してみます。
JuliaにはSet
と呼ばれる集合を表す型があります。これは型パラメータを取ることができる型で型パラメータをT
とすると、Set{T}
と書けます。つまりSet{T}
はT
型の要素を持つ集合です。このとき型パラメータT
はどのような型でも良いのですが(注: 実行時にはhash
を実装しているなどの制約がある)、逆にT
に何らかの制約(T<:Integer
など)を付けた型は作ることができませんでした。しかし、次期バージョンのJuliaの型システムでは、Set{T} where T <: Integer
などのようにT
型に制約をつけることができます。以下の例では、型制約がない場合とある場合でどのようにサブタイピングの結果が異なるかを示しています。
julia> Set{T} where T # 型制約なし Set{T} where T julia> Set{Int} <: (Set{T} where T) true julia> Set{Float64} <: (Set{T} where T) true julia> Set{T} where T <: Integer # 型制約あり Set{T} where T<:Integer julia> Int <: Integer true julia> Set{Int} <: (Set{T} where T <: Integer) true julia> Float64 <: Integer false julia> Set{Float64} <: (Set{T} where T <: Integer) false
サブタイピングに使えるということは、多重ディスパッチにも使えるということです。実際、以下の様な引数の型も書けます。
julia> f(xs::(Tuple{V,V} where V <: AbstractVector{T}), x::T) where T <: Integer = (xs[1] + x, xs[2] - x) f (generic function with 1 method) julia> f(([1,2,3], [1,2,3]), 2) ([3,4,5],[-1,0,1])
これが実際どのように役立つのかは、実例が無いので私にもまだよく分かっていません。また、記法に関してはまだ変更が入るかもしれません。
もうひとつ重要な新機能は、依存する関数の自動更新です。Julia v0.5では以下のように関数g
を実行した後f
を更新してもf
の変更は反映されませんでした。
julia> f(x) = x + 1 f (generic function with 1 method) julia> g(x) = f(2x) g (generic function with 1 method) julia> g(1) 3 julia> f(x) = x - 1 WARNING: Method definition f(Any) in module Main at REPL[1]:1 overwritten at REPL[4]:1. f (generic function with 1 method) julia> g(1) 3
これでは、開発中など頻繁に関数を書き換えるときにいちいち再起動する必要があります。Juliaはコードをコンパイルするためこのような動的な変更は実装不可能だと私は思っていましたが、現在進行中の変更(https://github.com/JuliaLang/julia/pull/17057)を取り入れると、f
の変更が正しくg
の実行に反映されるようになります。
julia> f(x) = x + 1 f (generic function with 1 method) julia> g(x) = f(2x) g (generic function with 1 method) julia> g(1) 3 julia> f(x) = x - 1 WARNING: Method definition f(Any) in module Main at REPL[1]:1 overwritten at REPL[4]:1. f (generic function with 2 methods) julia> g(1) 1
大変不思議です。どうやらworldという利用可能なメソッドをチェックできるようにする機構が組み込まれているようなのですが、私には詳しく分かりません。もし何か分かりましたら、また記事にしようと思います。
この記事ではJulia言語の現状を簡単にまとめてみました。Julia言語以外のコミュニティの動向などは、他の機会にまたご紹介しようと思います。明日は、Ishitaさんの担当です。
Julia言語の0.5の変更点
9月20日にJulia言語の最新版である0.5がリリースされました。Juliaのメーリングリストに投稿されたアナウンスメントはこちらです: https://groups.google.com/d/msg/julia-users/J2DiH1GnM8o/aO2Ku8o-CgAJ
きっと近いうちに本家のブログで詳しい変更点の紹介があると思いますが、私のブログでも一足先に主要な変更点をご紹介しようと思います。
クロージャの効率化
とりわけ重要な変更点として挙げられるのがJuliaのクロージャが効率化されたことです。
0.4までのJuliaでは、Juliaの関数に関数を渡したりJuliaで関数を返すような関数を作ると、その実効速度が極めて遅いことが問題でした。
これは、すべての関数が Function
という型にまとめられていたせいで、Juliaのコンパイラが特化したコードを吐けないせいでした。
これが0.5では関数は各々自分専用の型を持つように変更されたため、クロージャを使ったコードでも遅くなるということが無くなりました。
関数を受け取る関数の例としてよく上げられるのがmap
関数でしょう。
以下の簡単なベンチマークを見てもJulia 0.5では0.4に比べて10倍ほど速度が向上しています。
Julia 0.5.0:
julia> f = x -> 2x (::#1) (generic function with 1 method) julia> x = rand(100000); julia> map(f, x); julia> @time map(f, x); 0.000572 seconds (134 allocations: 789.078 KB)
Julia 0.4.6:
julia> f = x -> 2x (anonymous function) julia> x = rand(100000); julia> map(f, x); jjulia> @time map(f, x); 0.005628 seconds (200.01 k allocations: 3.815 MB)
クロージャを使うようなプログラムではこれでデザインの幅が広がり、今まで実行効率の観点からできなかったコードが書けるようになります。
ジェネレータ式
ジェネレータ式という新しい式が0.5から加わりました。
これは、イテレータをラップして式にしたようなもので、値を次々生成することができる式です。
簡単な例では、 2x for x in 1:9
のようなものが普通の値として扱えるようになりました。
julia> g = (2x for x in 1:9) Base.Generator{UnitRange{Int64},##5#6}(#5,1:9) julia> collect(g) 9-element Array{Int64,1}: 2 4 6 8 10 12 14 16 18
もちろん関数はこのジェネレータを値として受け取ることもできるため、次のような計算もできます。
julia> sum(g) 90
また、if
で値をフィルターすることも可能です。
julia> h = (2x for x in 1:9 if isodd(x)) Base.Generator{Filter{Base.#isodd,UnitRange{Int64}},##7#8}(#7,Filter{Base.#isod d,UnitRange{Int64}}(isodd,1:9)) julia> collect(h) 5-element Array{Int64,1}: 2 6 10 14 18 julia> sum(h) 50
文字列型の統合
Julia 0.4の標準ライブラリに大量にあった文字列型が整理され、0.5では主にString
型とSubString
型のみになりました。
String
型は従来のUTF8String
型に当たるもので、ユニコード文字列を表現できます。SubString
型は従来からありましたが、これはString
型の一部分を切り出す際に使われています。
将来的にはSubString
型もString
型に統合される予定のようですが、0.5時点では0.4のASCIIString
とUTF8String
がString
型に統合されたと見るのが良さそうです。
Julia 0.5.0:
julia> typeof("foo") String julia> typeof("いえい") String
Julia 0.4.6:
julia> typeof("foo") ASCIIString julia> typeof("いえい") UTF8String
昔の文字列型はLegacyStrings.jlパッケージに移動され、UTF-16などのエンコーディング方式はStringEncodings.jlパッケージで新たにサポートされるようです。
Fused broadcasting構文
Fused broadcasting構文は、Juliaのベクトル計算を効率化する新しい構文の拡張です。 0.4までのJuliaでは、配列に対して何度も関数適用を行うとその都度新しい配列が生成されていました。 新しいJuliaでは、これが構文レベルで融合(fuse)されるようになります。
sin.(cos.(x))
という式があったとしましょう。
0.4ではまずcos.(x)
が評価され、x
の各要素にcos
関数を適用した新しい配列が作られます。
続いてその各要素にsin
関数が適用され、また新しい配列が作られ、先ほどの配列は将来的にGCにより破棄されます。
0.5では、この式はまず broadcast(x -> sin(cos(x)), x)
のような式に変換されます。
これは、x
から値を一つずつ取り出してcos
とsin
を順に適用するため、一時的な配列の生成が起きず、一気に2つの関数を適用した配列が生成されます。
さらに、x .= ...
のような構文もbroadcast!(identity, x, ...)
に変換されるので配列への代入がin-placeに行うことができるようになりました。
マルチスレッドのサポート
マルチスレッドを使った計算の並列化が新たにサポートされました。
今のところ、@threads
マクロを使ってfor文の並列化を行うことができます。
例えば、以下のように複数のタスクを並列に処理することができます。
using Base.Threads function dotasks(tasks) @threads for i in 1:endof(tasks) dotask(tasks[i]) end end
配列などを領域で分割しておけば配列に対する処理の並列化もできますし、私の作ったBGZFStreams.jlではgzipの解凍処理をこのマルチスレッド機能を使って並列化しています。
テストフレームワークの強化
Julia標準のテストフレームワークがより強化され、テストをきれいに構造化できるようになりました。
今までBase.Test
から提供されていた機能では複数のテストをまとめる機能がなく、フラットにすべてのテストを書き連ねるしか方法がありませんでした。
しかし0.5からは@testset
マクロが加わり、次のようにテストを構造化することができるようになりました。
@testset "sum" begin @test sum([1,2,3]) == 6 @test sum([1.0, 2.0, 3.0]) == 6.0 end
その他の変更点
0.5では内部的なものも含め多数の変更がなされています。そのすべてはJuliaのリリースノートから参照することができますので、気になる方は是非調べてみて下さい。
JuliaからRを使う
先日のWACODE夏期講習でRCall.jlのデモをしたら、やはりウケが良かったようなので改めて紹介をします。
RCall.jlはJuliaからR言語の機能を呼び出すツールです。データの受け渡しからREPLでのインタラクティブな実行・プロットも簡単にできます。Juliaを使ってみたいけど、Rの豊富な資産を捨てる訳にはいかないといった方にはピッタリのライブラリです。
インストールは、Juliaの標準的な方法通り、julia -e 'Pkg.update(); Pkg.add(“RCall”)’
を実行して下さい。これで最新版のRCall.jlがインストールされることになります。尚、次期Juliaのリリース候補v0.5-RC1では現在動かないようですが、リリースブランチでは直っているのでRC2では使えると思います。
簡単な演算で正しくインストールできたかを確認しましょう。JuliaのREPLを起動して、using RCall
としてRCall.jlを読み込み、R”1 + 1”
と打ってR言語で1 + 1
を実行してみます。以下のように計算結果が表示されればOKです。
R”…"
という記法は、Juliaの非標準文字列と呼ばれる機能を利用したものです。この文字列の内側ではRのコードとして解釈され、Rの処理系がコードを実行してくれます。ダブルクウォートが含まれる場合はバックスラッシュでエスケープするか、R””” … “””
とトリプルクウォートを使うこともできます。
R側の評価値を取り出す場合には、rcopy
が使えます。
julia> rcopy(R"1 + 1") 2.0 julia> rcopy(R"1:5") 5-element Array{Int32,1}: 1 2 3 4 5 julia> rcopy(R"""list(x=10, y="foo")""") Dict{Symbol,Any} with 2 entries: :y => "foo" :x => 10.0
R側へJuliaの値を渡すには$
で変数を埋め込みます。
julia> x = 1 1 julia> R"1 + $x" RCall.RObject{RCall.RealSxp} [1] 2
JuliaのDataFrames.jlパッケージが提供しているDataFrame
なども自動的にRのdata.frame
へと変換してくれます。
julia> using RDatasets julia> iris = dataset("datasets", "iris"); julia> typeof(iris) DataFrames.DataFrame julia> R"head($iris)" RCall.RObject{RCall.VecSxp} SepalLength SepalWidth PetalLength PetalWidth Species 1 5.1 3.5 1.4 0.2 setosa 2 4.9 3.0 1.4 0.2 setosa 3 4.7 3.2 1.3 0.2 setosa 4 4.6 3.1 1.5 0.2 setosa 5 5.0 3.6 1.4 0.2 setosa 6 5.4 3.9 1.7 0.4 setosa
プロットも簡単にできます。以下のようにすると、Juliaから渡したデータをRがプロットしてくれます。
julia> R"""pairs($iris[1:4], pch=21, bg=c("red", "green3", "blue")[unclass($iris$Species)])""" RCall.RObject{RCall.NilSxp} NULL
いちいちR”…”
で囲んだりが面倒なら、RのREPLに入りましょう。JuliaのREPLから$キーを押すと、RのREPLに入れます。Juliaに戻るにはバックスペースです。
ここではRの文法が自由に使えるので、Rユーザーにとってはやりやすいでしょう。もちろん、Julia側からデータを受け取ることもR”…”
の時と同じように可能です。
R> library(ggplot2) R> ggplot($iris, aes(x=SepalLength, y=SepalWidth, color=Species)) + geom_point()
Rのセッションではヘルプなども参照できますし補完も効きます。Juliaは気になるけどRから離れるのはチョット…という方は是非一度試してみてください。