読者です 読者をやめる 読者になる 読者になる

プログラミングの実験場

Haskell/Webアプリ/画像処理/可視化/ITによる生産性向上 など

Julialangでの配列演算のベクトル化

任意の1変数関数、2変数関数はベクトル化ができ、行列の要素への適用が劇的に高速化する。

julia> A = rand(10000,10000);

julia> sq(x) = x*x;

julia> sq_vec(x) = x*x;

julia> @vectorize_1arg Float64 sq_vec
sq_vec (generic function with 4 methods)

julia> @time((B=sq(A);))
elapsed time: 26.948158237 seconds (816528604 bytes allocated, 0.06% gc time)

julia> @time((B=sq_vec(A);));
elapsed time: 0.520107978 seconds (800380984 bytes allocated)
  • 行末のセミコロンは、値の出力をしないという指示(Matlabと同様)。
  • @で始まるのはマクロ。@timeで式の実行時間を計測したり、@allocatedで式の実行時に確保したメモリサイズを表示したりできる。
  • @timeの中に入れている(B=sq(A);)というのはcompound expressionといい、複数の式をまとめることができる。(a=1; b=2; a*b)といった具合。

実行環境はMac mini 2012(Core i7 2.3 GHz)。要素数1億、800 MBのデータに対して0.5秒で済んだ。ベクトル化しない場合に比べて約50倍の高速化!

Julialangでの多次元配列使用のメモ: 配列の演算

基本の操作

添字による要素のアクセス

julia> a = [x*10+a for x=0:2, a=0:5]
3x6 Array{Int64,2}:
  0   1   2   3   4   5
 10  11  12  13  14  15
 20  21  22  23  24  25

julia> a[2,4]
13

julia> a[1:3,4]
3-element Array{Int64,1}:
  3
 13
 23

julia> a[2,2:4]
1x3 Array{Int64,2}:
 11  12  13

julia> a[2:4,4]
ERROR: BoundsError()
 in checkbounds at abstractarray.jl:65
 in getindex at multidimensional.jl:49

julia> a[2,:]
1x6 Array{Int64,2}:
 10  11  12  13  14  15

julia> a[2,3:end]
1x4 Array{Int64,2}:
 12  13  14  15

julia> a[2,[true,false,false,false,true,false]]
1x2 Array{Int64,2}:
 10  14

最後のBoolの配列による取り出しは、find()で変換したインデックスの配列を使っている。 julia> find([true,false,false,false,true,false]) 2-element Array{Int64,1}: 1 5

要素ごとの演算

  • Matlabのように、.+, .*などが使える。
julia> a = [x+a for x=0:2, a=0:5]
3x6 Array{Int64,2}:
 0  1  2  3  4  5
 1  2  3  4  5  6
 2  3  4  5  6  7

julia> a .% 2
3x6 Array{Int64,2}:
 0  1  0  1  0  1
 1  0  1  0  1  0
 0  1  0  1  0  1
  • 1変数のスカラー関数に行列を与えると自動で各要素への適用になる。
julia> f(x)=2x; f(a)
3x6 Array{Int64,2}:
 0  2  4   6   8  10
 2  4  6   8  10  12
 4  6  8  10  12  14

julia> g(x,y)=2x+y; g(a)
ERROR: `g` has no method matching g(::Array{Int64,2})
  • 2変数以上の関数はbroadcast(f,A, B, ...)という関数が使える。
julia> a = [x*10+a for x=0:2, a=0:5]
3x6 Array{Int64,2}:
  0   1   2   3   4   5
 10  11  12  13  14  15
 20  21  22  23  24  25

julia> b = [(x+y) % 2 for x=0:2,y=0:5]
3x6 Array{Int64,2}:
 0  1  0  1  0  1
 1  0  1  0  1  0
 0  1  0  1  0  1

julia> broadcast((x,y) -> y*x^(y+1), a, b)
3x6 Array{Int64,2}:
   0    1    0    9    0   25
 100    0  144    0  196    0
   0  441    0  529    0  625
  • In-placeでのbroadcast
julia> a = [x*10+a for x=0:2, a=0:5]
3x6 Array{Int64,2}:
  0   1   2   3   4   5
 10  11  12  13  14  15
 20  21  22  23  24  25

julia> broadcast! ((x,y) -> y*x^(y+1), a, a, b)
3x6 Array{Int64,2}:
   0    1    0    9    0   25
 100    0  144    0  196    0
   0  441    0  529    0  625

julia> a
3x6 Array{Int64,2}:
   0    1    0    9    0   25
 100    0  144    0  196    0
   0  441    0  529    0  625

Julialangでの多次元配列使用のメモ: 配列の作成

多次元配列に関するドキュメントなど

自分のための覚書、チートシート。Juliaのバージョンは0.3を使用。

公式ドキュメント

Multi-dimensional Arrays — Julia Language 0.3.6-pre documentation

これらも参照

http://sebastianraschka.com/Articles/2014_matrix_cheatsheet_table.html

100 julialang exercises - Qiita

http://bogumilkaminski.pl/files/julia_express.pdf

初期化

zeros, ones, eye

julia> zeros(3)
3-element Array{Float64,1}:
 0.0
 0.0
 0.0

julia> zeros(3,1)
3x1 Array{Float64,2}:
 0.0
 0.0
 0.0
 
julia> zeros(3,4)
3x4 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> ones(4,2)
4x2 Array{Float64,2}:
 1.0  1.0
 1.0  1.0
 1.0  1.0
 1.0  1.0

julia> eye(3)
3x3 Array{Float64,2}:
 1.0  0.0  0.0
 0.0  1.0  0.0
 0.0  0.0  1.0

julia> eye(Int,3)
3x3 Array{Int64,2}:
 1  0  0
 0  1  0
 0  0  1

julia> eye(3,4)
3x4 Array{Float64,2}:
 1.0  0.0  0.0  0.0
 0.0  1.0  0.0  0.0
 0.0  0.0  1.0  0.0

julia> eye(3,4,2)
ERROR: `eye` has no method matching eye(::Int64, ::Int64, ::Int64)
  • 2次元のNx1 Arrayと1次元のN-element Arrayは別の型を持つ。
  • 最初の引数で数値の型を指定できる。デフォルトはFloat64
  • eyeは正方行列じゃなくても使える。

rand

julia> rand()
0.8480348007358423

julia> rand(1,3)
1x3 Array{Float64,2}:
 0.327055  0.251911  0.37398

julia> rand(1:4,1,10)
1x10 Array{Int64,2}:
 3  4  3  4  1  3  2  3  2  3

julia> rand(1.0:0.3:9.0,1,10)
1x10 Array{Float64,2}:
 1.9  5.5  1.0  1.3  2.2  5.2  6.4  2.2  4.9  2.5
  • 最初の引数で型あるいは乱数の範囲を指定できる。 範囲は以下の様な指定方法がある。
1:4 # == 1, 2, 3, 4
1:2:10 # == 1, 3, 5, 7, 9
1:1.5:10 # == 1, 2,5, ...8.5, 10.0

範囲型の値は[]で囲むことで1次元の配列に変換できる。

julia> [1:1.5:10]'
1x7 Array{Float64,2}:
 1.0  2.5  4.0  5.5  7.0  8.5  10.0

'は転置行列を作る後置演算子。

内包表記

  • 以下の式において、x,a,zがそれぞれ指定された範囲を動き、x+a+zを各セルの値とする。
  • 変数名はなんでも良い。変数の登場順(この場合、x, a, z)に行、列、第3次元となる。
julia> a = [x+a+z for x=0:2, a=0:5, z=0:1]
3x6x2 Array{Int64,3}:
[:, :, 1] =
 0  1  2  3  4  5
 1  2  3  4  5  6
 2  3  4  5  6  7

[:, :, 2] =
 1  2  3  4  5  6
 2  3  4  5  6  7
 3  4  5  6  7  8

Julialangでの画像処理

Juliaの概要

Juliaという科学技術計算言語が最近出てきた(2009年に開発開始、2012年にオープンソース)のだが、これがかなり良さそうに見える。MatlabやRがやっていることを置き換えられるようなものであると思う。言語仕様やインフラがモダンという印象で、プログラマにとっての使い勝手も良いためか急速に人気も出てきているようで、現時点で中央レポジトリに500程度のパッケージができている。大学の講義などでも使われ始めている。

言語仕様と処理系として特徴的なのは

  • 実行速度が高速。このベンチマークを見ると分かるように、Javaより速く、Cの2倍程度に収まる。Juliaでループを書いてもMatlabのように遅くならない。
    • Juliaの処理系自体(Githubのレポ)も6割以上Juliaで書かれている。
  • ノイズが少なく読みやすい文法。文法に変な制限なども無くRubyのような感じで自然に書ける。Unicodeも使用可能。
  • 多重ディスパッチによる多相。
  • 値や関数の引数のオプショナルな型アノテーション(高速化にも寄与する)。
    • なお、変数には型はない。違う型の値を再代入してもエラーにならない。constはある。
  • レコード型によるオブジェクトの作成。
  • オブジェクトを作れる具象型と、作れない抽象型があり、具象型間の継承は無い。具象型は複数の抽象型を継承しうる。
  • 型パラメータによるジェネリクス。型パラメータには数値も使え、実行時にチェックされる。行列の次元を型で区別したりできる。
  • オンラインでブラウザ上で動くREPL(対話環境)やIDEもある。
  • Githubを使ったパッケージ管理システム。自作のパッケージの登録も簡単にできるようだ。
  • 外部とのバインディングが簡単。
    • 外部ライブラリの利用がバインディングがCで書かれた外部のコードの呼び出しもバインディング不要でccallという関数一つで呼べる。
    • Juliaの関数を外から呼び出すのも簡単そうに見える。
  • マクロ・メタプログラミング機能がある。
    • プログラム中で型情報も含めた構文木(S式)を作成・操作し、実行可能。

多重ディスパッチというのは聞いたことがなかったが、関数の実行時に全引数の型に応じて異なる実装が呼ばれる仕組みで、関数はオブジェクトに所属するのではなく、トップレベルの存在になる。

メタプログラミング機能も、凝ったことができる自由度が無限大という意味ではなかなか興味をそそられる機能だと思う。

関連リンク

現状Pythonで科学技術計算している人も多いと思うが、Pythonとの比較を含めた詳しい解説は、こちらを参照。

PythonistaのためのJulia100問100答 - りんごがでている

速度など、Pythonとの比較。

Julia vs Python: ビットコインオプションのモンテカルロシミュレーション | once upon a time,跡地

高い抽象度とネイティブレベルの壁を超えるという好意的な評価。

「私がJuliaを推す理由」(翻訳) - アドファイブ日記

具体例

私がやりたいことの実例をコードを挙げて例示してみる。

  • 3チャンネルの画像(別々のファイル)をROIで切り取った後、RGB画像として保存する。
import Images: Image,RGB,Gray,imread
import FixedPointNumbers: UfixedBase

function main(channel :: String = "561multi")
    files = split(readall("image_list.txt"),'\n')
    rois = readdlm("roi_list.tsv", '\t')
    idx = findfirst(a -> a == "Slice", rois[1,:])
    count = 0
    for (f,i) in zip(files,[1:length(files)])
        if chomp(f) == ""
            continue
        end
        roilist = choose(rois[2:end,:], idx, i)
        println("$(size(roilist,1)) rois.")

        img1 = imread(replace(f,"CH","488nm"))
        img2 = imread(replace(f,"CH","561nm"))
        img3 = imread(replace(f,"CH","640nm"))

        for i=1:size(roilist,1)
            count += 1
            print(".")
            ns = roilist[i,:]
            crop_with_roi(img1, img2, img3, ns, rois[1,:], count)
        end
    end
end


function choose(dt::Array{Any,2}, idx::Int, slice::Int)
    res = dt[1,:] :: Array{Any,2}
    for i=2:size(dt)[1]
        if dt[i,idx] == slice
            res = [res; dt[i,:]]
        end
    end
    res
end

typealias GrayImage Array{Gray{UfixedBase{Uint16,16}},2}

macro debug(v)
    return :(println($v); $v)
end

function autocontrast(im::GrayImage)
    mx = maximum(im)
    mn = minimum(im)
    (im .- mn) ./ (mx-mn+0.001)
end

function mk_color(img1::GrayImage,img2::GrayImage,img3::GrayImage)
    a,b = size(img1)
    cr2 = Array(UfixedBase{Uint16,16},a,b,3)
    cr2[:,:,1] = autocontrast(img1)
    cr2[:,:,2] = autocontrast(img2)
    cr2[:,:,3] = autocontrast(img3)
    convert(Image{RGB},cr2)
end

function crop_with_roi(img1,img2,img3,row,header,num)
    out_folder = "/Volumes/MacintoshHD/output_test/"
    y = convert(Int,row[findfirst(a -> a == "BY", header)][1])
    x = convert(Int,row[findfirst(a -> a == "BX", header)][1])
    h = convert(Int,row[findfirst(a -> a == "Height", header)][1])
    w = convert(Int,row[findfirst(a -> a == "Width", header)][1])
    cr1 = img1[x+1:x+w,y+1:y+h]
    cr2 = img2[x+1:x+w,y+1:y+h]
    cr3 = img3[x+1:x+w,y+1:y+h]
    Images.imwrite(mk_color(cr1,cr2,cr3), out_folder * "test$num.tif")
end

CSVの読み方がよくわからなかったので適当なコードを書いた。Imagesというパッケージを追加で使っている。JuliaのREPLから、Pkg.add("Images")で入れられる。

多次元配列の使い方はMatlabのような感じ。ただし、いくつか違いもある。

Noteworthy Differences from other Languages — Julia Language 0.4.0-dev documentation

250x250ピクセルの3チャネルTIFFを300個切り出して保存するのに20秒ほど。二回目以降は事前コンパイルされているせいか、7秒程度に短縮された。悪くない。

若干不満な点

Juliaにはかなり期待していてこれから書く画像処理はしばらくJuliaを使ってみようと思うのだが、若干不満な部分もあった。

関数の型がない

関数はすべてFunctionという一つの型であり(例えば、typeof(map) == Function)、Haskellみたいにmap :: (a -> b) -> [a] -> [b]みたいな型アノテーションは書けない。 異なる型の関数を同じ名前に割り当てる多重ディスパッチとは相性が悪いようにも思うので、しょうがない気もするが。

実行時エラーがそれなりに多い

型が強すぎる言語だと、特に書き捨ての小さなスクリプトでは面倒になるので、要はバランスなのだが、演算中に型が合わずにエラーということが感覚としては結構多かった。慣れの問題はあると思うが。 変数の未定義なども実行時エラーとなるので、そのあたり実行前にもうちょっと強くチェックしてくれるとありがたいと思う。

importがかなり遅い

Pkg.add("Images")でパッケージを入れた上に、importするのだが、REPLを起動した後の最初のimportが異様に遅い。数秒待たされる。何かやり方が間違っているのかも...。

Scala.jsでD3.jsを使う

Scala.jsの概要

Scala.js http://www.scala-js.org/ という、ScalaからJavaScriptのコードを生成できるシステムがあり、ちょっといじって見ているが、かなり良く出来ているという印象を持った。IntelliJの補完も使って、型安全+Scalaの抽象化能力を最大限に活用できる。

http://tototoshi.github.io/slides/sendagaya-js-scala-js/#1

このスライドで書かれているように、初出は去年なのだが、進歩が猛烈なスピード。現在ver. 0.5.1。この調子で行けば、じきに高機能AltJSの代表格となるのではないだろうか。

使ってみて感じたメリット。

  • Scalaの言語仕様のほぼ全機能(リフレクション以外)を使える。
    • 型安全で、Scalaのパワフルな抽象化が使える。JSのオブジェクトも型付けできる。
    • 当然IntelliJで補完も効く。
  • JavaScriptとの接続(Scala.js←→JSの双方向)が明快、便利。
    • これまで使った感触だと、あまりハマリポイントはない模様。
    • ソースマップも生成される。
  • 最大最適化後のファイルサイズは150KB程度まで小さくなった(最小のプロジェクトで)。
    • 手持ちのMac Mini (Late 2012, Core i7 2.3GHz)でビルド時間は10秒ほど。
    • 最適化を弱めると1MBのJSファイルが出来た。ビルド時間は3秒。
  • 導入が簡単。
    • Scalaプロジェクトの.sbtファイルに2行書き足すだけですぐに使える。
      • スタンドアロンのコンパイラも最近用意された。
    • HTMLからの呼び出しは、生成された単一のJSファイルを読み込むだけ。
    • DOM操作やjQueryのバインディングも用意されている。

Scala.jsの導入

@mizchiさんの記事(http://mizchi.hatenablog.com/entry/2014/02/14/211400)でも書かれているように、まず短いサンプルとしてサンプルのレポジトリをクローンして走らせてみるのが良い。

git clone --depth 1 https://github.com/lihaoyi/workbench-example-app
cd workbench-example-app
sbt gen-idea
sbt ~fastOptJS

これでIntelliJでスラスラ書きながら、変更を検出して自動コンパイルできる。

チュートリアルは短くて分かりやすい。 http://www.scala-js.org/doc/tutorial.html

言語はScalaそのものなので、Scalaの言語以外に追加で必要な知識はJSとのFFI(言語間接続)だけで、それさえわかればすぐに使える。

FFIのポイント
Scala.js→外部のJS
  • JSのオブジェクトはjs.Objectトレイトを継承したトレイトとしてScala上にマッピングされる。
外部のJS→Scala.js
  • js.Array(2,3,4,5)というようにしてJSネイティブの配列を作れる。
  • js.Dynamic.literal("name" -> "John", "age" -> 19)というようにしてJSネイティブのオブジェクトを作れる。
  • js.FunctionNはscala.FunctionNと等価。コールバックに用いる。

Scala.js→D3.jsのFFI

型付きのFFIを作成したいわけだが、TypeScript用に書かれたD3.jsの型定義があるので、それを元に以下のコンバータを使ってScala.js用の定義を作成してやる。このレポジトリには他にも多数のライブラリがある。

D3.d.ts: https://github.com/DefinitelyTyped/DefinitelyTyped/tree/master/d3

Importer from TypeScript type definitions to Scala.js: https://github.com/sjrd/scala-js-ts-importer

コンバータが吐き出したScalaファイルは大方ちゃんと出来ているようにみえるが、実際走らせようとすると、コンパイルエラーと実行時エラーがあり、いくつか修正する必要があった。 以下に生成されたもの、手動修正後のソース、そしてその2つのdiffを示しておく。すべてを手で修正する時間と気力はなかったので、とりあえずサンプルコードが動くようになるところまでの最低限の修正をした(後ほど完全なバージョンを作って公開したい)。

Importerが生成したソース: https://gist.github.com/nebuta/7beb3d6e8d17a7ba5142

修正後のソース: https://gist.github.com/nebuta/c20e0db85bf3e7a0dab1

Diff: https://gist.github.com/nebuta/c703bc0317822c253752

修正のポイント:

  • overrideを追加する必要があるメソッドがいくつかある。(IntelliJ上で表示されるので分かる)
  • D3のSelectionオブジェクトはJSの配列を拡張しているので、js.Arrayをextendsする。
  • js.Any型の引数のコールバックを使う関数(例えばattr()など)は、以下のようにsubtypeも含めたシグネチャに変える必要がある。
def attr[A <: Any](name: js.String, valueFunction: scala.Function2[A, js.Number, Any]): Selection = ???
  • js.Array[A]は(どういうわけか)invariantで定義されている、配列を取る関数はdef func[A <: HogeType](x: A)としてやる。そうしないと、型が合わないことが頻発する。
  • Scalaの仕様上、オーバーロードしている関数で、複数の関数にデフォルト変数を設定することができない。適宜= ???を削除する。

使用例

作成したバインディングを用いて、D3.jsの公式サイトにある例をScala.jsに書き換えたものを以下に挙げる。

デモ: http://nebuta.github.io/d3js-scala/

ソースコード: https://github.com/nebuta/ScalaJS-D3-Example/tree/master/src/main/scala/example

まとめ・感想

  • ウェブアプリのサーバーとクライアントを同じ言語で書きたいというのはわりと根強くある希望だと思うが、例えばScalaでPlay frameworkをサーバーに使い、Scala.jsをクライアントに使う、というふうにすれば、かなり生産性が高くなりそう。
  • AltJSとJSとの相互運用は、双方の変数・関数がどのようにマッピングされるかというセマンティクスの接続(?、そういう言葉があるか知らないが)を理解するのが概して難しい。
    • AltJSの中だけでなるべく閉じたような構造にすれば制御しやすくなる。
    • 違う言語を接続するのは本質的にチャレンジングな課題。CoffeeScriptやTypeScriptのような「薄い」ラッパーであれば比較的単純になるが。
    • 最近の大型JSフレームワーク(AngularJS, Ember, Meteorなど)との組み合わせも、どのくらいスムーズに統合できるのか、興味がある。
  • エラーを取り除くために、バインディングのコードにアドホックともいえる修正を結構したので、まだ安心感がない...。(FFI一般にそうだが、正しく書いていないと不可解な実行時エラーが生じうる。)

余談だが、D3.jsについての感想:

  • D3.jsは公式サイトに載っている実例の鮮やかさも手伝ってか、人気が非常に高い。しかし、いくつか改良点も考えられる。
    • コンセプトを理解するのがわりと難しい。
    • メソッドチェインが長く、何が操作対象なのかが理解しづらい時がある。
      • このことが、型付き言語でD3を使いたくなる動機である。
      • 以前作った、HaskellからD3のコードを生成するライブラリ。複雑なJSライブラリへの型付けの実験でもある。http://hackage.haskell.org/package/d3js
    • D3.js自体が抽象化レイヤーであるためか、他の抽象化と相性が悪いような気がする。
    • JSフレームワークとD3.jsを組み合わせようと思いいろいろ調べたのだが、両方ともデータ/ビューバインディングを実現するライブラリであるため、統合がスムーズにできない感が強い。
      • 典型的には、JSフレームワークによる画面更新が終わった時に、データバインドをオフにしたようなDOM要素にD3を使って描画する、という継ぎ接ぎ的な構造になる。

ということで、強力なAltJSが出てきている今、そろそろD3.jsに影響された、型付きのAltJSを使ってより高度に構造化した次世代可視化ライブラリが出てきても良い頃かと思われる。Scala.jsやPureScript、あるいはHasteを使って一から作るのもいいかもしれない。

参考資料

他にも役に立ちそうなスライド http://www.slideshare.net/kamekoopa/kyonkekkonlt

D3.jsの基本パターン

D3.jsをいじっているが、最初に期待していたのよりも難しい。

公式サイトには資料が結構多くあるのだが、あまり整理されていない感じがある。

サイトに載っている実装例が大変に華麗でmotivatingなのだが、その分ちょっと複雑なソースコードが多い。

なので、以下で基本操作を最短でマスターするための道のりについての覚え書き。

データとDOM要素の対応付け

まず、データとDOM要素(主にSVG)を対応させる方法を理解する必要がある。以下のパターンがデータからDOM生成する基本である。

var dat = [{r: 10, x: 100, name: "Alice"},{r: 8, x: 150, name: "Bob"},{r: 15, x: 200, name: "Chris"}];

var sel = d3.select("body").append("svg").attr('width',300).attr('height',300);

sel.selectAll("circle")
    .data(dat, function(d){return d.name;})
    .enter()
    .append("circle")
    .attr('fill','pink')
    .attr('r',function(d){return d.r;})
    .attr('cx', function(d){return d.x;})
    .attr('cy', function(d,i){return i*50;});

第2文はjQueryみたいな感じで直感的である。selはbody内に追加されたsvg要素を示す。

attr()の第2引数には値として、定数あるいはコールバック関数を指定できる。 コールバック関数は

  • 第1引数: データポイント(数値、オブジェクト、配列)の配列
  • 第2引数: 配列のインデックス

を受け取る。

第3文のシーケンスはちょっとあれ?となると思う。これを理解するには、以下の記事が役に立つ。

selectAllというのは仮想的なプレースホルダをselの中に生成する操作である。selectとselectAllの違いはここにある。一見ちょっとややこしいが、selectAllは階層構造を潰して平たい配列を作るが、selectは階層構造を保持する、という違いがある。

enter(), exit()という関数は、data()でデータを割り当てたあとにメソッドチェインして、DOMを更新する関数。

ここで、データには、data()の第2引数でコールバック関数を渡すことで、identity keyを設定することができる。このkeyが対応するDOM要素にも記憶されるので、データをdata()で再度割り当てた時に、DOM要素をそれに応じて更新できる。

もし割り当てた新しいデータが、選択された既存のDOM要素と一致しない場合:

  • enter().append('DOM要素') を呼ぶと、足りない分のDOMが追加される。
  • exit().remove() を呼ぶと、新しいデータに対応物が存在しないDOM要素が削除される。

このデータとDOM要素のマッチングの際に、上記のidentity keyが使われる。もしdata()でidentity keyを指定していなかった場合は、単に配列のインデックスが使われる。つまり、data(dat,function(d,i){return i;})としたのと同じことになる。

それを例示しているのが、以下の例

トランジション

transition()をメソッドチェインの中に挟むと、それ以降のattr()やstyle()での変更がアニメーションになる。また、

  • delay()でアニメーション開始時刻を設定できる。
  • duration()でアニメーションにかかる時間を設定できる

この2つとも引数として関数を指定できるので、データポイントに依存した設定が可能。以下に例を示す。

   sel.transition()
    .delay(function(d,i){return 1000*i;})
    .duration(1000)
    .attr('r',function(d,i){return d.r*2});

スケーリング

データの値と描画座標の間の変換を設定できる。これを有効活用すれば、自分でデータから描画座標を計算する必要はあまり無いはず。 線形、指数、対数のスケールや、座標の離散化、quantileでの座標変換を設定することもできる。 また、数値→座標だけでなく、

  • 時間から座標を計算するtime scale や
  • 系列データから座標を計算するordinal scale

というのもある。

その他

  • http://bl.ocks.org/mbostock/4062085

    • nest()を使ったデータをグループ化する方法が分かる。
  • グラフなどを書くための高レベルAPI

    • d3.layout.***の関数が用意されている。

Halide事始め

MIT Halide(http://halide-lang.org/)というC++の高速画像処理ライブラリが去年リリースされた。C++内のDSLを使っていて、画像処理を、

  • 計算のアルゴリズム
    • アルゴリズム自体は副作用がなく、純粋関数的に定義される。
  • 実装の詳細(「スケジュール」)
    • スケジュールは、計算の並列化と、複数の計算ステップの融合を最適化する方法を指定する。

を分離した形で簡潔に書ける。

凄いのは、抽象的に書いたアルゴリズムからライブラリによって低レベルコードを出力させたほうが、手書きでチューニングしたコードよりも速いことすらあること(という作者の主張)。たとえば最近開発されたlocal Laplacian filterというの複雑なアルゴリズム(http://people.csail.mit.edu/sparis/publi/2011/siggraph/、pixelwise operationやstencil operationだけでなく、任意の画素にアクセスするような計算も含む。)が、OpenMPとIntel IPPで並列化した手書きのC++よりも2.1倍速く、コードが3.7倍短い(http://people.csail.mit.edu/jrk/halide12/)とか。

このHalide、かなり有望な感じがするのだけど、なんせ資料が少ない。公式サイトにもサンプルコードがあるのみで、ドキュメントらしいドキュメントはない(追記:http://halide-lang.org/Halide/index.html にAPIレファレンスがあった)。しかしこの技術、DSL、画像処理、並列計算、純粋関数型の計算による抽象化、と僕の興味に見事にマッチしている。そこで以下に覚書を兼ねてライブラリの概要を記す。

論文

使い方

Githubのレポhttps://github.com/halide/Halideをクローンしてmakeでビルドする。llvm (>= 3.2)とllvm-configが使えるようにしておく必要がある。(私はここで微妙にハマった...。)

ごく簡単なチュートリアルがtutorialsフォルダの中にあるが、情報量は少ないので、むしろtestとappsの2つのフォルダ内のコードが参考になりそう。

  • 実行時にJITコンパイルする方法
  • C++を実行時にDSLを展開した別のC++を生成する方法

の2種類があるが、使える機能は同じようだ。単に試す分には簡単に使えるJITコンパイルで良い。

Halide::FuncとHalide::Varの2つのクラスを使って、例えば以下のように画像を表現できる。ソースを読んでいないので仕組みはよくわからないがこれはなんかマジックっぽい感じがする。

    Halide::Func output;
    Halide::Var x, y;
    output(x, y) = x + y;
    output.realize(800,600);

tutorialsフォルダにあるサンプルコードをコンパイルするには、以下の様な感じで良い。

$ clang++ -lHalide `libpng-config --cflags --ldflags` lesson_02.cpp

Funcクラスのインスタンスとして定義した関数に対してrealize()を呼ぶことで計算を実行する。 並列化などの「スケジュール」を指定する場合は、realize()の前にそれを指定する。

testフォルダ内にたくさん小さなコードがあるが、簡単な例は以下のようなものがある。

  • two_vector_args.cpp
  • vector_extern.cpp
  • vector_bounds_inference.cpp
  • lambda.cpp
    • ラムダ式が使えるので、その例。

これまでの例はpixelwiseの操作だが、線形フィルタなど、近傍のピクセルを使った計算(stencil computation)もできる。

実際に使ってみる

以下の様なコードをサンプルコードからコピペ+適当に変更して書いた。画像の明るさを調整するプログラム。 https://gist.github.com/nebuta/6149137

コンパイルして実行してみる。

$ clang++ -I ../include -L ../bin -lHalide `libpng-config --cflags --ldflags` brighter.cpp
$ ./a.out 
Input: 1536 x 2560 x 3
Start measurement.
Vectorized vs scalar (average over 1000 times): 2.20 ms 23.02 ms. Speedup = 10.482
Success!

約4メガピクセルの画像が2msで処理できている。良い感じだ。スケジュールをたった1行指定するだけで10倍に加速した。ちなみにスケジュールの指定を適切にしないと速くならないし、何もしないよりも遅くなることすらあった。

追記:画像の入出力、pointwise operation(コントラスト調整・2値化など), stencil operation(ガウスフィルタなど)といった基本操作をサンプルコードからコピペ+適当に変更してまとめてみた。https://gist.github.com/nebuta/6152664 コンパイルにlibpngが必要。 実行結果は

$ clang++ -lHalide `libpng-config --cflags --ldflags` halide-examples.cpp 
$ ./a.out 
1.337 ms (1000 times average): pointwise
0.104 ms (1000 times average): pointwise,parallel
12.9x speedup.
8.153 ms (10 times average): image_rgb_to_gray
4.168 ms (10 times average): image_rgb_to_gray, par
2.0x speedup.
56.128 ms (10 times average): blur
9.090 ms (10 times average): blur,par
6.2x speedup.

良い感じだ。

まとめ

MIT Halideを使うメリットは、

  • 画像処理のアルゴリズムが純粋関数的に抽象的・簡潔に書くことが可能で、実装の詳細(「スケジュール」)は、アルゴリズムから分離して後から最適化できる。
  • スケジュールによって、単一ステップの並列化と複数ステップの融合を同時に行うことができる。
  • 簡単なスケジュールの指定で、劇的に速度が向上する。

現状不足している点は、

  • スケジュールを人間が考えて手書きしないとならず、それだけちょっと難しい気がする。
  • ドキュメントがほとんどなく、サンプルコードを読んでライブラリの機能を推測する必要がある。
  • 特に、スケジュールの書き方に関する指針がほとんど示されていない。

Halideのようなアプローチが一般化すれば、OpenCVのような既存の大きなフレームワークに頼ることなく、自分の必要な機能を柔軟かつ簡潔にボトムアップで構築することも十分可能になると思われる。

2報目の論文ではスケジュールの自動最適化も行なっているので、そのうちユーザー向けにも実装されることが期待できる。複雑なアルゴリズムの場合、スケジュールの組み合わせを単純にランダムに探索すると組み合わせ爆発してしまうそうで、遺伝的アルゴリズムで効率的な探索をしている。