新しいブログを作りましたのでお知らせします

新しいブログはhttps://www.bicycle1885.org/です。 こちらは今後更新しない予定です。 よかったら最初の記事をご覧ください。 www.bicycle1885.org

Juliaのfor文はなぜ速い?

1年以上ぶりのブログ更新ですね。

atsuoishimotoさんのブログ記事で，私のブログ記事が言及されていました。 atsuoishimoto.hatenablog.com

書いてある内容はまさにその通りで，for文を回すこと自体が取り立てて遅いわけではなく，その中で何らかの処理をする際に動的型付け言語であるPythonの性質上避けて通れない処理があり，その処理をfor文の中で何度も繰り返すことで結果としてfor文全体の処理が遅くなるという内容でした。

そうすると，ひとつの疑問として，なぜ同じ動的型付け言語であるはずのJulia言語では同様の処理が高速なのかが気になるところでしょう。この記事ではそのあたりを少し解説しようと思います。

Juliaの型推論

問題となっているコードは，以下のように２つのベクトルの内積を計算するコードです。

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

Python:

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通りの方法があるよ!

1. 標準ライブラリのreadcsv
2. DataFrames.jlのreadtable
3. 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でやろうとしている人達が出てきて、いくつかパッケージができているよ。

• DataFramesMeta.jl
• Query.jl
• StructuredQueries.jl

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から離れるのはチョット…という方は是非一度試してみてください。