2014-07-06

JuliaTokyo #1 を開催しました

Julia

日本で(多分)最初のJuliaの勉強会 JuliaTokyo #1 を開催しました。イベントページはこちらです。 http://juliatokyo.connpass.com/event/6891/

当日は40人ほどの参加者が集まり、いつの間にJuliaがそんなに一般化したんだという感じでしたが、Julia歴1週間とか当日の朝から始めましたというように、何処からかJuliaが熱いらしいとの噂を聞きつけてJuliaを使ってみようという新しいもの好きの人たちが多い印象でした。

私はパッケージの作り方についてお話をしてきた次第です。

以下は、それぞれの発表と私のひとこと感想です。

メインセッション

Juliaのこれまでとこれから - @sorami

Juliaの背景的な情報やどんなエコシステムになっているかを紹介していました。 JuliaはPythonやRからユーザを取り込みつつ、共に進化していく道を辿るようです。

Julia100本ノック - @chezou

Julia 100 exercises #JuliaTokyo from chezou

Juliaでnumpyの100 numpy exercisesを書いた話です。 Juliaのやり方とnumpyのやり方はベクトル化などの使い方がかなり違い、苦労があったようです。あとJuliaのベクトル内包表記は速いのでガンガン使っていきましょう！

Juliaのパッケージを作ろう! - @bicycle1885

Juliaのパッケージをつくろう! from Kenta Sato

私の発表です。 Juliaのパッケージをどのようにつくるかを解説しました。パッケージ名やディレクトリ構成などについて、ドキュメントに書いてない慣習などが大事だったりします。

ビジュアライズJulia - @nezuq

ビジュアライズJulia from nezuQ

"The Grammar of Graphics"のlayerやscaleの概念や、Gadfly.jlとの関係の話でした。可視化関係のライブラリはまだそれほどなので、今後充実してくると良いですね。

メカ女子将棋 - @kimrin

メカ女子将棋Julia tokyo#1 from Takeshi Kimura

Juliaで将棋AIを作った話でした。やっぱりJuliaはCなどに比べ開発がずっとしやすいですね。盤面の表現にはJuliaの128bit整数を使っているようです。

A Brief JuliaCon Report / Natural Language Processing with Julia - Pontus Stenetorpさん

https://github.com/JuliaCon/presentations/blob/master/JuliaNLP/JuliaNLP.pdf

先月Chicagoで開催されたJuliaConのレポートと、そのときに発表されたJuliaでの自然言語処理のお話をしていただきました。 JuliaConでは、アカデミックの人がメインだったものの、実際にサービスにJuliaを投入している会社もあったようです。 Juliaで自然言語のパーサallenを作って、パフォーマンスも満足いくものになったようです。

LT

Juliaでちゃんとしたコードを書く - 関さん

当日の朝からJuliaを始めて、どうするとJuliaが遅くなるかのお話でした。マニュアルのperformance tipsマジ重要です。読みましょう!!

Julia0.3でランダムフォレスト - @gepuro

Julia0.3でランダムフォレスト from Atsushi Hayakawa

DecisionTree.jlを使ってみたようです。 DecisionTree.jlはID3を決定木構築のアルゴリズムに使っているので、CARTが使いたければ私のRandomForests.jlを使いましょう！

Plotly Julia API - Yasunobu Igarashiさん

Plotly Julia API from Yasunobu Igarashi

Plotlyというデータの可視化や共有のためのプラットフォームの紹介でした。 JuliaのAPIもあるんですね、すごい！

LightTableを使いましょう - @QuantixResearch

LightTable + Jewel でJuliaを開発しようという話でした。インラインで変数の値が見れたり、補完もかなり強力なようです。あとで使ってみよう！

感想

新しいおもちゃを手にしたときの子供のような、Juliaに対する期待と興奮が垣間見れたように思います。おぉ！Juliaすごい！いいじゃん！というような熱気がありました。 Twitterでしか見たことない人と実際にお話したり、情報交換をしたりも大変有益でした。反省点としては、ちょっと初心者向けの話が少なすぎたかなと思います。次回は初心者向けセッションを設けて欲しいとの要望が多かったです。会場を用意していただいた株式会社ブレインパッドの皆様、本当にありがとうございました。

あと、本家Juliaステッカーをもらえたのがすごい嬉しかった！

本場シカゴからやって来たステッカー貰った(((o(*ﾟ▽ﾟ*)o))) pic.twitter.com/1AJeQCp0ml
— 佐藤ヾ(⌒(｢ 'ω')｢ｶﾞｵｰ建太 (@bicycle1885) 2014, 7月 5

2014-06-10

HaskellのFusionがあれば速度と抽象化を両立できる

Haskell

データの分析をする際には配列やベクトルは欠かせないデータ構造です。大体どの言語にも大体配列は用意されていて、そこにサンプルのデータ等を入れて統計量を計算したり関数に渡して回帰をしたりするわけです。ベクトルという単位でデータの塊を扱うものの、実際のアルゴリズムではベクトル内の要素一つひとつを見ていって何か処理をしたり計算をすることが多いでしょう。その際、命令型言語ではforループを陽に使って要素にアクセスすることになります。

簡単な例を見てみましょう。ベクトルv1とv2の内積を求める関数をC言語で書くと以下のようになります。

function dot(double* v1, double* v2, size_t n)
{
    double s = 0.0;
    for (int i = 0; i < n; i++) {
        s += v1[i] * v2[i];
    }
    return s;
}

特に難しいところはありませんね:)

ベクトル化された記法

しかし、大きなベクトルに対してこうして陽にループを書くことは、PythonやRではある種の「悪手」であると認識されています。理由は単純で、PythonやRのインタープリタのループが遅すぎるため、NumPyのようなCやFortranなどで書かれた高速なライブラリに配列(実際にはPythonのオブジェクトが持つメモリ上のバッファへのポインタ)を渡して計算させるのが定石になっています*1。

NumPyやSciPyではこうしたベクトルに対する関数や演算子が大量に用意されており、これらを組み合わせれば、明示的なforループを避けて簡潔に高速な計算ができるようになっているわけです。こうした操作をベクトル化(vectorized)された操作と言われています。例えば、先程と同様２つのベクトルの内積を計算するならNumPyにnumpy.dot関数が予め用意されています。

import bumpy as np

v1 = np.array([1.0, 2.0, 3.0])
v2 = np.array([2.0, 3.0, 4.0])

print(np.dot(v1, v2))

ベクトルの要素ごとの積(*)と総和を求める関数を使えば、内積計算を以下のようにも書くことができます。

def dot(v1, v2):
    return np.sum(v1 * v2)

分散の計算

こうしたベクトル化された関数を組み合わせて標本分散(不偏分散)を計算してみましょう。数式での定義は以下のとおりです。

f:id:bicycle1885:20140610232956p:plain

これを単純にNumPyを使ったPythonのコードに写してみましょう。(もちろん、NumPyには分散を計算するnp.var関数が用意されています！)

def var(x):
    return np.sum((x - np.mean(x))**2) / (np.size(x) - 1)

ベクトル化された関数を使えば、定義式をほとんどそのまま移すだけで、分散を計算する関数が作れました。

しかし、この方法には問題があります。それは、計算の途中に出てくる式の計算のために不要な一時配列が確保され、余計なメモリを食っていることです。上の例では、x - np.mean(x)では、元の配列xと同じ大きさの配列が確保され、平均からの差分をそこに格納していく形になっています。(x - np.mean(x))**2ではさらに二乗の結果を格納する配列が確保されています。これらの一時的な配列は、本来分散を計算するためにはまったく不要なものです。内積の例でも、np.sum(v1 * v2)のところでnp.sumに渡す一時ベクトルが新たに作られてしまっています。

def var(x):
    return np.sum((x - np.mean(x))**2) / (np.size(x) - 1)
#                  ~~~~~~~~~~~~~~
#                  ^
#                 ~~~~~~~~~~~~~~~~~~~
#                 ^
#                 余計なベクトルが2個作られている

Pythonのforループは遅くて使えないため、これを回避するためには、CやCythonで陽にループを書くハメになります。

#include <stdlib.h>
#include <stdio.h>

// Two-pass algorithm:
//   http://en.wikipedia.org/wiki/Algorithms_for_calculating_variance#Two-pass_algorithm
double var(double* v, size_t n)
{
    double mean = 0.0;

    for (int i = 0; i < n; i++) {
        mean += v[i];
    }

    mean /= n;

    double ssd = 0.0;

    for (int i = 0; i < n; i++) {
        ssd += (v[i] - mean) * (v[i] - mean);
    }

    return ssd / (n - 1);
}

先ほどのベクトル化された計算と比べると、ループに展開したコードはかなり読みにくく、定義式との対応がわかりづらくなってしまいました。

ちなみに、Juliaでは、マクロを使ってベクトル化された記法からループへと展開するライブラリDevectorize.jlがあります。

fusionを使って宣言的に書き、効率的なコードを吐く

宣言的なベクトル化された記法をしつつ、ループと同等の効率的な実行を可能にするのがfusionです。 fusionでは、中間にできてしまう一時的なベクトルをライブラリによる自動的なコードの書き換え規則とコンパイラによる最適化で排除し、効率的なコードを生成する技術です。 Haskellのvectorパッケージでは、ほとんど意識する必要なくベクトル計算のfusionを行ってくれます。 unboxedな値を格納するベクトルを使うなら、Data.Vector.Unboxedモジュールをインポートして使います。ここでは見た目を簡潔にするため、Preludeの名前が衝突する関数を隠しています。

import Prelude hiding　(sum, zipWith)
import Data.Vector.Unboxed

dot :: Vector Double -> Vector Double -> Double
dot v1 v2 = sum $ zipWith (*) v1 v2

v1 * v2ではなくzipWithを使ってる点が若干分かりにくいかもしれませんが、ループを書かずに宣言的に内積が定義できています。それでいて、実はこのdot関数は計算の中間データを保持するベクトルを作っていません!

fusionがどうやって動いているのか

実際にfusionがどのような仕組みになっているか簡単な例で見てみましょう。ここでは、Stream Fusion - From Lists to Streams to Nothing at Allのサンプルの一部ををGHCで試せるように書いてみます。これは、stream fusionというfusionの一例です。 vectorパッケージではVector型に対してfusionを行っていますが、ここではリストを使っています。

minimal stream fusion

例えば、map f . map gというコードは関数合成の結合則とstream . unstream ≡ idという関係式から、

map f . map g ≡ (unstream . mapS f . stream) . (unstream . mapS g . stream)
              ≡ unstream . mapS f . (stream . unstream) . mapS g . stream
              ≡ ustream . mapS f . mapS g . stream

となります。ここで、streamとunstreamが打ち消しあうことで、リストとStreamの変換が消えますが、これはGHCのrewrite ruleによって実現しています。さらに、GHCの最適化によりStepの生成も排除され、効率的なコードが吐かれるわけです。リストに対するstream fusionのより詳しい実装は論文の著者等によるstream-fusionパッケージを参照してください。ただ、最近のGHCは既にリストののfusionをサポートしているため、上のコードやstream-fusionパッケージを使ったものは普通に書いたものと同等のパフォーマンスになりました*2。

分散の計算、再び

分散の計算をvectorパッケージを使って書いてみましょう。 fusionの効果を見るため、わざと関数を分けて定義しています。

import Prelude hiding (sum, length, zipWith, map)
import System.Environment (getArgs)
import Data.Vector.Unboxed as V
 
-- | pairwise product
(.*) :: Vector Double -> Vector Double -> Vector Double
x .* y = zipWith (*) x y
 
-- | sample mean
mean :: Vector Double -> Double
mean v = sum v / fromIntegral (length v)
 
-- | deviations from the mean
deviation :: Vector Double -> Vector Double
deviation v = map (\x -> x - mean v) v
 
-- | vectorized square
square :: Vector Double -> Vector Double
square v = v .* v
 
-- | sum of squared deviations from the mean
ssd :: Vector Double -> Double
ssd = sum . square . deviation
 
-- | unbiased sample variance
var :: Vector Double -> Double
var v = ssd v / fromIntegral (length v - 1)

このようにかなり宣言的にプログラムを書くことができ、それぞれの関数は単独で使うことも、他の場所で部品として使うこともできるなど再利用性も極めて高いです。Cの方だと可能なのはせいぜい平均値を計算するmean関数をくくりだすくらいでしょうか。

パフォーマンス比較

大きい配列に対してパフォーマンスを比較しみましょう。長さ100,000,000の配列に対してCのループを使った実装とパフォーマンスを比べてみます。使用したコードは以下のリンクにあります。 Gist - variance

Haskell: 0.49s user 0.35s system 99% cpu 0.839 total
C:       0.48s user 0.33s system 98% cpu 0.815 total

Haskellの実装とCの実装とでほとんど同じ速度がでていますね。Haskellのコンパイルには-fllvmでLLVMのコードを吐くようにすると何割か速くなりました。

本当に一時データが割り当てられていないかも見てみましょう。 Doubleは8byteなので、長さ100,000,000の入力ベクトルに対して800,000,000バイト程度が必要になります。以下のヒープに割り当てられた領域の大きさからも、不要な配列が割り当てられていないのが分かります。

% ./variance-hs $(cat n) +RTS -s -RTS

     800,121,472 bytes allocated in the heap
          10,104 bytes copied during GC
          44,312 bytes maximum residency (1 sample(s))
          21,224 bytes maximum slop
             764 MB total memory in use (0 MB lost due to fragmentation)

                                    Tot time (elapsed)  Avg pause  Max pause
  Gen  0         1 colls,     0 par    0.00s    0.00s     0.0001s    0.0001s
  Gen  1         1 colls,     0 par    0.00s    0.06s     0.0562s    0.0562s

  INIT    time    0.00s  (  0.00s elapsed)
  MUT     time    0.50s  (  0.72s elapsed)
  GC      time    0.00s  (  0.06s elapsed)
  EXIT    time    0.00s  (  0.06s elapsed)
  Total   time    0.50s  (  0.83s elapsed)

  %GC     time       0.0%  (6.8% elapsed)

  Alloc rate    1,610,816,342 bytes per MUT second

  Productivity 100.0% of total user, 60.0% of total elapsed

新しいfusion

このようにfusionは強力なツールですが、今までのstream fusionにはいくつか限界があります。最後に、こうした制約を取っ払うために発案された2つのfusionを紹介しましょう。実装の詳細に関しては、私の能力の限界を超えるため、興味のある方は元の論文を参照してください。

Generalized Stream Fusion

2つのベクトルを結合して新しいベクトルを作りたいことはよくあります。しかし、既存のstream fusionでは一度に1つづつしか値を取り出せないため、ベクトルを塊として扱うmemcpyのような効率的なコピーができません。また、同じ理由で複数の値に一気に操作をするSIMD命令も使うことができません。

しかしGeneralized Stream Fusionでは素朴なStreamをより柔軟なBundleという単位にまとめることでこれらの問題を解決してくれます。 2つのベクトルの内積を計算するベンチマークでは、SSEを使うようにしたCのコードと同等かそれ以上の速度を出しています*3。

f:id:bicycle1885:20140610233645j:plain

Generalized Stream Fusionでは、既存のHaskellのコードをほとんど変更することなく、合成可能な再利用性を保ったまま高速化できるのが魅力です。この機能はvectorパッケージ上に実装されています。現在OpenBLASやEigen3などの数値計算ライブラリと同等のパフォーマンスは出せていませんが、開発が進めばこれらのライブラリに迫るものができるかもしれません。

Data Flow Fusion

もう一つのstream fusionの拡張が、Data Flow Fusionです。元のstream fusionではベクトルの値を変換して別のベクトルなり集計値を計算する「消費者」は常にひとりだけに限られています。しかし、1つのベクトルから複数の統計量を計算したりする分岐をしたいことが常です。大きなベクトルを複数回ループすることはキャッシュの有効利用ができないため、1つのループにまとめるのに比べて不利になってしまいます。

以下のコードは入力ベクトル(vec1)の要素をすべてインクリメントし(vec2)、正数のみを集め(vec3)たものとその中の最大の数をタプルにして返す関数です。

filterMax :: Vector Int -> (Vector Int, Int)
filterMax vec1 =
  let vec2 = map    (+ 1) vec1
      vec3 = filter (> 0) vec2
      n    = foldl max 0 vec3
  in (vec3, n)

この処理は一度のループで実行可能なはずですが、stream fusionでは1つのfusionにできない処理のようです。また、複数のベクトルに対して同時にループを行うときにループカウンタが重複してしまいレジスタを無駄遣いしてしまうこともパフォーマンスを悪化させてしまいます。

こうした問題をベクトル長の変化やデータフローを解析する高度なfusionによって解決するのがData Flow Fusionです*4。実装は、Repaの配列に対するGHCの最適化プラグイン(repa-plugin)として提供されているようです。この最適化により通常のstream fusion(Stream)に比べてdata flow fusion(Flow)では大きく性能が向上し、人間が手でfusionしたCのコード(Hand-fused C)に迫るパフォーマンスの向上を実現しています。

f:id:bicycle1885:20140610233649j:plain

fusionは既に実用的な技術になっており、より強力になったfusionが使えるようになるのもそう遠くないかもしれません。

*1:Pythonで数値計算のコツ：for文書いたら負けかなと思っている

*2:GHC User’s Guide - 7.21.4. List fusion

*3:Exploiting Vector Instructions with Generalized Stream Fusion

*4:Data Flow Fusion with Series Expressions in Haskell

2014-06-08

Pythonのコマンドライン引数処理の決定版 docopt (あとJuliaに移植したよ)

Python Julia

Pythonをよく使う人にはよく知ってる人も多いのですが、docoptという便利ライブラリがあります。 docoptはargparseやoptparseのようなコマンドライン引数をパースするライブラリなのですが、その発想がコロンブスの卵なのです。

例えばPython標準のargparseだと、argparseのAPIを組み合わせてパーサを組み立てるわけです。するとパーサと一緒にヘルプも作ってくれて、"program --help"などとすると自動生成されたヘルプを表示してくれるようになります。しかし、そのAPIを覚えるのが大変で、毎回ドキュメントを読まないと忘れちゃうわけです。

import argparse

parser = argparse.ArgumentParser(description='Process some integers.')
parser.add_argument('integers', metavar='N', type=int, nargs='+',
                   help='an integer for the accumulator')
parser.add_argument('--sum', dest='accumulate', action='store_const',
                   const=sum, default=max,
                   help='sum the integers (default: find the max)')

args = parser.parse_args()
print(args.accumulate(args.integers))

https://docs.python.org/dev/library/argparse.html

その問題をdocoptではスマートに解決してくれます。docoptでは、パーサからヘルプを生成するのではなくて、ヘルプからパーサを生成します。そしてコマンドライン引数は辞書(dict)を使って取得できるようになります。 docoptモジュールがエクスポートしているのはdocopt関数のみで、複雑なAPIを覚える必要もなく、馴染みのあるヘルプを文法に従って書けばパーサを生成してくれます。

先の例をdocoptで書き換えれば、以下のようになります。

"""Process some integers.

usage: prog.py [-h] [--sum] N [N...]

options:
    -h, --help  show this help message and exit
    --sum        sum the integers (default: find the max)
"""

from docopt import docopt

if __name__ == '__main__':
    args = docopt(__doc__)
    Ns = map(int, args["N"])
    if args["--sum"]:
        print(sum(Ns))
    else:
        print(max(Ns))

とても直感的ではありませんか！使用すると、以下のようになります。引数が文法に合わなければ usage が表示されます。

/Users/kenta/tmp% python3 prog.py 4
4
/Users/kenta/tmp% python3 prog.py 4 3 5 2
5
/Users/kenta/tmp% python3 prog.py --sum 4 3 5 2
14
/Users/kenta/tmp% python3 prog.py
Usage: prog.py [-h] [--sum] N [N...]
/Users/kenta/tmp% python3 prog.py -h
Process some integers.

usage: prog.py [-h] [--sum] N [N...]

options:
    -h, --help  show this help message and exit
    --sum        sum the integers (default: find the max)

ヘルプの文法は標準的なもので、以下の様なものがあります。

要素	説明	補足
`<argument>`	位置引数	`<x> <y>`に引数`1 2`を与えれば、`{"<x>": 1, "<y>": 2}`と結び付けられる
`--option`, `-o`	オプション	ハイフンで始まる引数
`(arguments)`	必須の引数	引数はデフォルトで必須だが、グループ化するときに使える
`[arguments]`	必須でない引数	`(...)`と違って、なくてもエラーにはならない
`arg1 \| arg2`	arg1もしくはarg2	3つ以上連ねることもできる (例. `--up \| --down \| --left \| --right`)
`args...`	可変数引数	返り値はリストになる

これらを組み合わせて、以下のように複雑なものも書けます。これをargparseで書くのは至難の業でしょう。

Naval Fate.

Usage:
  naval_fate ship new <name>...
  naval_fate ship <name> move <x> <y> [--speed=<kn>]
  naval_fate ship shoot <x> <y>
  naval_fate mine (set|remove) <x> <y> [--moored|--drifting]
  naval_fate -h | --help
  naval_fate --version

Options:
  -h --help     Show this screen.
  --version     Show version.
  --speed=<kn>  Speed in knots [default: 10].
  --moored      Moored (anchored) mine.
  --drifting    Drifting mine.

(from docopt.org)

docoptを使うモチベーションやより詳しい文法については公式サイト(http://docopt.org/)や作者によるPyCon UK 2012の動画を御覧ください。

PyCon UK 2012: Create beautiful command-line ...

言語非依存なのも利点の1つです。一度使い方を覚えれば、他言語でもAPIを特別覚える必要がなく使えます。オリジナルのPythonの他にも、Haskell, Go, C, Ruby, CoffeeScriptなどにも移植されています。

Julia版のdocopt - DocOpt.jl

個人的に最近Juliaをよく使っているのですが、docoptが無かったので移植したしました。 GitHubに上げたのは何ヶ月か前なのですが、最近やっと正式なdocoptファミリーの一員になり、Juliaの正式なパッケージにしたので改めて紹介します。

ソースコードはこちらです。https://github.com/docopt/DocOpt.jl

文法や使い方はPython版と同じですので、既にPython版を使ったことのある方は何も覚えることはありません。 using DocOptしてdocopt関数をインポートし、ヘルプを渡すだけです。

doc = """Naval Fate.

Usage:
  naval_fate.py ship new <name>...
  naval_fate.py ship <name> move <x> <y> [--speed=<kn>]
  naval_fate.py ship shoot <x> <y>
  naval_fate.py mine (set|remove) <x> <y> [--moored|--drifting]
  naval_fate.py -h | --help
  naval_fate.py --version

Options:
  -h --help     Show this screen.
  --version     Show version.
  --speed=<kn>  Speed in knots [default: 10].
  --moored      Moored (anchored) mine.
  --drifting    Drifting mine.

"""

using DocOpt  # import docopt function

arguments = docopt(doc, version=v"2.0.0")
dump(arguments)

パース結果はPython同様、辞書になります。

インストールはJuliaの対話環境から、Pkg.add("DocOpt")とすればインストールできます。残念ながら今のところJuliaのv0.2系をサポートしていないので、v0.3のコンパイル済みprerelease版をインストールするか、ソースからビルドすることになります。 v0.2で使いたいという要望があれば、v0.2もサポートしようと思います。

2014-04-11

Haskellでデータ処理がしたい！

Haskell

この記事はIIJの@kazu_yamamotoさんの下で、アルバイトとして書いています。最終的な目的はHaskellでの統計処理やデータ処理のツールを整備しようというものですが、その下調べをした内容をこの記事にまとめています。ご意見などコメントいただけるとありがたいです。

今から3年前の2011年のICFPに、Emily G. Mitchell 氏のHaskellを使ったシミュレーションに関するExperience Reportが掲載されました。

[PDF] Functional Programming through Deep Time - Modeling the ﬁrst complex ecosystems on Earth

このレポートでは筆者が行ったエディアカラ紀の生態系シミュレーションの実装に関してプログラミング言語の選択や、使った言語にどのような利点・欠点があったかが述べられています。初めはRを用いていたようですが、実行速度とシミュレーションのコード修正をするのに問題があってHaskellを選択した経緯とそのときに感じたHaskellでデータ処理をする際の問題点はとても興味深いものでした。筆者は元々Haskellの使い手だったわけではないようですが、夫があのNeil Michell氏であることもあってか、シュミレーションにはHaskellを使い、プロットなどを含むデータ処理にはRを使ったようです。

その中でRと比較したHaskellの利点として以下のようなことが挙げられています。

構文が一貫していること
純粋さと静的型付けのおかげでコードのリファクタリングが容易であったこと
型システムのおかげで、デバッグが容易であったこと

しかし、Haskellまわりのツールには不満点もあるようで、

異なるパッケージの似たようなデータ型で互いにデータを変換するのが困難であること
Haskellパッケージのインストールが難しいこと (RはGUIでパッケージをインストールすることもできる)
標準ライブラリで統計量の計算やプロットができないこと

などが挙げられています。

現状Haskellの統計パッケージあたりがどういうのがあるか調べてみると、Haskell Wikiにあるのは、

の3つだけです。(hstatsは2008年で止まっている)

プロットに関しては、

といった感じです。正直、充実しているとは言い難い状況ですね。

現状の機能では大変なことを実感するには、簡単な具体例を考えてみれば分かりやすいでしょう。以下の例は、2012年のHaskell cafeにあったスレッドの一部を翻訳したものです。

練習として、GHCiを使って適当な浮動小数点数の表になっているCSVファイルを読み込んで、ある列と別の列の間で線形回帰をして、回帰直線と一緒に散布図をプロットしてみてください。残差の正規性もチェックしてみるといいかも。大きなデータも扱えるようにするには、bytestringとかattoparsecなんかも使う必要があるでしょう。x86_64のGHCiからhmatrix/gslの関数を使うと、セグメンテーションフォルトとかバスエラーになる既知のバグがあることにも気をつけましょう*1。何か当たり前なことを見落としてるかもしれませんが、(データの)コンテナ・永続化・パース・統計・プロットをするのにどのパッケージを選ぶのにすっっごく時間がかかるでしょう。

[Haskell-cafe] Mathematics and Statistics libraries

確かにこれをHaskellでやるのはちょっと考えてみただけどもやりたくない感じですね。まずCSVのパースなのでcassavaあたりを使うといいのでしょうか。cassavaのData.CSVモジュールを見るとdecodeというのがあるのでこれを使えば良さそうですね。でも引数がByteString型なのでbytestringパッケージのData.ByteStringモジュールをインポートしてファイルを読む必要がありますね。返り値はVector型なのでそれから必要な列を取り出して線形回帰すればいいのですが、Haskellで線形回帰をするライブラリを探してぐぐってみるとstatistics-linregがありますね。よし、これをcabalでインストールしてとかやってられません。 しかしRやPythonなら、慣れている人なら5分とかからずにやってしまうでしょう。試しにRでやってみるとたった6行です。

library(ggplot2)  # ggplot2ライブラリを読み込む
data = read.csv("data.csv")  # ヘッダーつきCSVファイルを読み込む
data.lm = lm(y ~ x, data)  # 線形回帰 
plot(data.lm)  # 残差の正規性をQ-Qプロットでチェック
ggplot(data, aes(x=x, y=y)) + geom_point() + stat_smooth(method="lm")  # 試しにプロット
ggsave("data.png")  # プロットを画像ファイルにセーブ

入力のCSVファイル:

id,x,y
0,0.125,2.140
1,0.065.567
2,-0.079,3.115
3,-0.193,1.765
...

プロット:

f:id:bicycle1885:20140411183421p:plain

この簡潔さを可能にしているRの特徴として挙げられるのは、

表データを扱う組み込みのデータ構造 (データフレーム)
データ処理のための組み込み関数 (read.csv(), lm() など)
プロットのためのライブラリ (ggplot2, latticeなど)

などでしょう。 Pythonでも組み込みの機能はRほどではないにしろ、データフレームにはpandas、データ処理にはnumpy/scipy/scikit-learn、プロットにはmatplotlibがあり、IPythonを使えばRと同様インタラクティブに処理可能になります。

これを真似て処理するデータを収めるコンテナを統一的に扱えるようにし、データ処理によく使う関数をひとつのモジュールをインポートするだけで使えるようにすれば、結構Haskellでもいけるようになるのではないでしょうか。加えて、快適に使うために欲しいものとしては、ドキュメントがサクサク読めて、データを眺めたりプロットをはめ込めるRStudioやIPython notebookのような統合的対話環境があると嬉しいです。

これらの機能を整備したとして、Haskellを使う利点は何でしょうか。先の論文の要旨も踏まえて私の考えを3点述べると、

静的な型付けによる安全性とメンテナンス性の向上
データ並列による処理の高速化
入出力の抽象化とストリーミング処理

これらはRやPythonではそう簡単には実現できない特徴ですので、うまくデータ処理のニーズと合致すればとても魅力的なツールになるはずです。最後に、データ解析者やHaskellerの方々の目線から、このような機能は必須だとか、Haskellのこのライブラリを使うといいなどの意見をコメントいただけるととてもありがたいです。

*1:GHCの7.8.xでは修正されているようです。 https://twitter.com/kazu_yamamoto/status/454741146053255168

2014-03-16

Juliaのデストラクター

Julia

Juliaのマニュアルを読むとコンストラクターの記述はちゃんとあるのですが、デストラクターに関しての記述が見当たりません。ググっても情報が出てこないようなので、ここにその方法を残しておきます。

JuliaにはC++の~ClassName()のようにデストラクター専用の関数はありませんが、コンストラクターの中でオブジェクトがGCされたときの処理を記述できます。やり方は簡単で、Base.finalizer()関数がまさにオブジェクトが破壊された時の処理を登録するための関数ですので、これを使うわけです。

finalizer(x, function)はxが到達不可能になったときにfunction(x)を呼ぶ仕組みになっているため、functionに資源を解放するような処理を置いておけば良いことになります。

何かDBなどとの接続を管理する型Connectionの例で考えてみると、以下のように何かmallocなどの直後にfinalizer()で後処理を登録すればいいことになります。

type Connection
    conn::Ptr{Void}
    open::Bool

    function Connection()
        this = new(c_malloc(sizeof_connection_info), false)
        finalizer(this.conn, connection_finalizer)
        err = ccall(
            (:connection_init, "libsomedb"),
            Cint,
            (Ptr{Void},),
            this.conn
        )

        if err != 0
            error("failed to create a new connection")
        end
        this.open = true
        this
    en

Juliaのソースで実際に使われている例を見てみましょう。正規表現のライブラリで確保したメモリへのポインタをfinalizer()で開放しています。

type Regex
    pattern::ByteString
    options::Uint32
    regex::Ptr{Void}
    extra::Ptr{Void}

    function Regex(pattern::String, options::Integer)
        pattern = bytestring(pattern)
        options = uint32(options)
        if (options & ~PCRE.OPTIONS_MASK) != 0
            error("invalid regex options: $options")
        end
        re = compile(new(pattern, options, C_NULL, C_NULL))
        finalizer(re,
            function(re::Regex)
                re.extra != C_NULL && PCRE.free_study(re.extra)
                re.regex != C_NULL && PCRE.free(re.regex)
            end)
        re
    end
end

https://github.com/JuliaLang/julia/blob/0f54ca52e350c89d821069cbdc349f1e49e43cbc/base/regex.jl#L7-27

注意する点としては、

finalizer()はimmutableな型には適用できない
コンストラクターでは最後に自分自身を返すことを忘れない

ことが挙げられます。

参考

2014-02-14

matplotlib入門

Python

matplotlibはPythonでグラフを描画するときなどに使われる標準的なライブラリです。画像ファイルを作るばかりでなく、簡単なアニメーションやインタラクティブなグラフを作ることも可能です。実際の例はmatplotlibサイトのギャラリーで見ることができます。

matplotlib/gallery

matplotlibは本家のサイトやどこかのブログにあるチュートリアルや例を描画してみるぶんには簡単なのですが、実際に自分でプロットするとなると基礎的な概念を理解していないと使いにくいライブラリでもあります。また、基礎的な概念を理解していないとドキュメントを参照する際にもどこを見て、どう実用すればいいのかわかりません。そこで、この記事ではそのあたりのmatplotlibの基礎を解説していきます。なお、Python自体の知識はある程度仮定していますが、matplotlib自体の実装に関わるようなことまでは解説しませんので Pythonの関数やクラス、モジュールなどがわかれば十分でしょう。描画するデータの作成にはnumpyを利用しています。

なお、この記事はまだ書いてる途中です。今後もっと網羅的かつ手短にmatplotlibを把握できる記事にしていくつもりです。

2012-12-24

Haskell Tips (文字列編)

Haskell

この記事はHaskell Advent 2012のための記事です。(遅れてすいません(ヽ´ω`))

Haskellで実際にプログラムを書く上では、様々なパッケージにある型や関数を利用するのが不可避になります。そういった便利ツールのうち文字列まわりについて調べたところを、Haskell最近始めたって人に紹介するのがこの記事の目的です。Haskell詳しい方には新しい発見は無いかもしれません。ごめんなさい。

この記事では、主に文字列に焦点を当てていますが、そのうちまた別のテーマでも書こうと思っています。

なお、前提としてるバージョンはHaskell Platform 2012.4.0.0と、それに付随するGHCやパッケージです。 OSも、Mac OS XやLinuxなどUNIX系の環境を想定しています。Windowsは持ってませんのでわかりません。ごめんなさい。

文字列

Haskellでは文字列の表現も色々あって、Haskellを始めたばかりの人はきっと混乱するでしょう。自分も大混乱しましたよ。

Haskellで文字列を扱うのによく使われる型は大きく分けて以下の３つです。

String
ByteString
- Strict
- Lazy
Text
- Strict
- Lazy

ByteString型はbytestringパッケージで、Textはtextパッケージから提供されてます。どちらもHaskell Platformに同封されてます。

ひとまず確認して行きましょう。

String

あんまり説明は要らないと思いますが、Stringはただの文字型Charのリストです。なのでlengthやmapなどのリスト操作関数はそのまんま使えて便利ですし、自然な発想です。

Stringでは一文字で５ワードも占めるので、メモリを喰いますし、速度もイマイチです。

でもStringはマルチバイト文字を扱えます。以下のようなものはランタイムのロケールを考慮して、正しく出力されるはずです。

putStrLn "鬼"

そうそう、入出力のエンコーディングなどのロケール情報は、GHC.IO.Encodingモジュールの getLocaleEncodingやgetLocaleFileSystemEncodingで取得できます。

import GHC.IO.Encoding

sample1 :: IO ()
sample1 = print =<< getLocaleEncoding

出力

UTF-8

他にもGHC.IO.Encodingにはエンコーディングを設定する機能などがあるようなので、この機能も必要があるかもしれません。

ByteString

「すごいHaskell」や「Read World Haskell」などを読まれた方はおなじみかと思います。bytestringパッケージで提供されているByteString型と関数で、効率良く文字列（というよりバイト列）を扱うことができます。

ByteStringが実際に収めているデータはWord8(8bitの符号なし整数)型なので、それより大きい要素を収めようとすると問題があります。マルチバイト文字など8bitに収まらないデータがそうです。

ByteStringには、実は同名の２つの型があります。ひとつはData.ByteString.Internalモジュールで定義されている正格評価のByteStringで、もうひとつはData.ByteString.Lazy.Internalモジュールにある遅延評価のByteStringです。内部的には、Lazyな方のByteStringは、StrictのByteStringをリストにしているだけです。

-- LazyなByteStringの定義
import qualified Data.ByteString.Internal as S
data ByteString = Empty | Chunk {-# UNPACK #-} !S.ByteString ByteString

正格と遅延の型を字面だけ見て組み合わせようとしても、型エラーが出るので気をつけましょう(｀･ω･´) 外部のライブラリを使うときは、ByteStringとだけ書いてあってもどちらか分からないので、ソースを読むなり、HTML版のHaddockのリンク先を確認するなりしないといけませんね。

じゃあ、これらの型の相互変換はどうするのかというと、bytestringパッケージのバージョン0.10.0.1では、 toStrictとfromStrictが、LazyとStrictを相互変換するための関数が定義されています。それ以前のバージョンでは、これらの関数は自作することになります。このへんの話は以下のリンクを参考にして下さい。

他に、Data.ByteStringにはモジュール名の末尾にChar8がつくData.ByteString.Char8とData.ByteString.Lazy.Char8もありますが、これらはByteStringが内部でもつWord8のデータとCharの相互変換をしてくれる関数を間にかませているだけなので、本質的には何も変わりません。変換の分のオーバーヘッドはありそうですが。

このへんの挙動をちょっと見てみましょう。

下のコードの頭にあるOverloadedStringsは、文字列リテラルをIsString型クラスのインスタンスへと自動的に変換してくれます。ですので、文字列リテラルに関しては、bytestringやtextパッケージで用意されているpack関数を明示的に呼ぶ必要はありません。

コメントでは文字に対応する文字符号化形式のコードを書いています。

{-# LANGUAGE OverloadedStrings #-}

import qualified Data.ByteString.Char8 as B8

sample :: IO ()
sample = do
    -- '鬼':
    --   UTF-16: 9B3C
    --   UTF-8:  E9ACBC
    -- '<':
    --   ASCII: 3C
    B8.putStrLn (B8.singleton '鬼')

    -- '平':
    --   UTF-16: 5E73
    --   UTF-8:  E5B9B3
    -- 's':
    --   ASCII: 73
    B8.putStrLn (B8.singleton '平')

    B8.putStrLn (B8.pack "鬼平")
    B8.putStrLn "鬼平"      -- OverloadedStringsがなければ、この書き方は型エラー

出力

<
s
<s
<s

漢字などマルチバイト文字は、Char型としての文字リテラルではUTF-16でエンコードされている様子です。しかし、ByteStringは内部的にはWord8(8bit)単位でしか扱えないので、 ByteStringに変換する際に、CharのデータはWord8(8bit)に切り詰められてるのが分かると思います。

文字列リテラルで指定されたマルチバイト文字を無理やりByteStringで扱おうとすると、特にエラーも出さず間違った結果を返すので注意が必要ですね。なので、ByteStringは文字列というよりは、生のバイト列を収めるデータ型と思ったほうが良いと思います。

IOを介すなどして外部から与えられたマルチバイト文字は、エンコーディングなど考慮されずにそのままバイト列として扱えます。なので、以下のようにユーザーの入力を受け取る場合は、入力がマルチバイト文字を含んでいても問題ありません。入力をそのまま吐いてるだけです。

sample :: IO ()
sample = do
    putStrLn "What is your name?"
    putStr "> "
    name <- B8.getLine
    B8.putStrLn $ "Hello, " `B8.append` name
    putStrLn $ (show $ B8.length name) ++ "-byte"

出力

What is your name?
> 鬼平
Hello, 鬼平
6-byte

バイト数がUTF-8でエンコードしたときにあたる6バイトになっちゃってますね。これでは文字数は正しく数えられませんし、文字境界も分からないのでgrepのようなことも出来ません。

Text

じゃぁこういったマルチバイト文字をバイト列でなく、ちゃんと文字列として扱うにはどうすればいいのかと言うと、textというパッケージがありまして、こいつを使うというわけです。

ちゃんと文字列として扱うとは、文字列の長さが正しく数えられるとか、特定の文字を探せるとかそういう話です。

textパッケージでは、文字列操作のためのモジュールとIOのためのモジュールが別のモジュール(Data.Text.IO)にあります。下のように、同じTという別名をつけてimportすると楽ちんです。

{-# LANGUAGE OverloadedStrings #-}

import qualified Data.Text    as T
import qualified Data.Text.IO as T

sample :: IO ()
sample = do
    T.putStrLn (T.singleton '鬼')
    T.putStrLn (T.singleton '平')

    T.putStrLn (T.pack "鬼平")
    T.putStrLn "鬼平"       -- これもOverloadedStringsが必要

出力

鬼
平
鬼平
鬼平

良い感じですね(´ω｀) Data.Text.IOではロケールを考慮してくれています。

先ほどのByteStringのときと同じサンプルも、Textなら大丈夫です。

sample :: IO ()
sample = do
    putStrLn "What is your name?"
    putStr "> "
    name <- T.getLine
    T.putStrLn $ "Hello, " `T.append` name
    putStrLn $ (show $ T.length name) ++ "-letter"

出力

What is your name?
> 鬼平
Hello, 鬼平
2-letter

ちゃんと文字数を数えています。バイト数じゃなくって。でも、Textのlength関数の計算量はO(n)です(´・ω・｀)

Data.Text.IOには他にファイルやHandleとText型のデータをやり取りする関数も用意されてます。

TextはByteStringとの相互変換も可能ですし、Lazyなバージョンもあります。相互変換に必要な関数はData.Text.Encodingモジュールにあります。

利用できるエンコーディングは以下の３つです。

UTF-8
UTF-16
UTF-32

下の例では、マルチバイト文字をByteStringで表現して、decode系の関数がそのByteStringをバリバリ喰ってTextにガガッと変換してくれます。

ちなみに、任意のバイト列を作成したい時は、以下の様に文字列リテラルの中でバックスラッシュ()と16進数を表すxで指定することができます。 10進数では、xは要りません。

{-# LANGUAGE OverloadedStrings #-}

import qualified Data.ByteString          as B
import qualified Data.Text                as T
import qualified Data.Text.IO             as T
import qualified Data.Text.Encoding       as T
import qualified Data.Text.Encoding.Error as T

sample :: IO ()
sample = do
    -- '鬼':
    --   UTF-8:  E9ACBC
    --   UTF-16: 9B3C
    let oniUTF8    = "\xE9\xAC\xBC" :: B.ByteString
        oniUTF16BE = "\x9B\x3C"     :: B.ByteString -- ビッグエンディアン
        oniUTF16LE = "\x3C\x9B"     :: B.ByteString -- リトリエンディアン

    T.putStrLn $ T.decodeUtf8    oniUTF8     -- Output: 鬼
    T.putStrLn $ T.decodeUtf16BE oniUTF16BE  -- Output: 鬼
    T.putStrLn $ T.decodeUtf16LE oniUTF16LE  -- Output: 鬼

    -- エンディアンの間違い!!
    T.putStrLn $ T.decodeUtf16LE oniUTF16BE  -- Output: 㲛

UTF-8ではバイトオーダーの指定は不要なのですが、 UTF-16では必要で、バイトオーダーを間違えると別の文字として認識されていしまいます。外部から受け取ったUTF-16でエンコードされているはずのByteStringなどを、Text型に変換する際はバイトオーダーに気をつけましょう。

あと、外部から受け取ったデータなどをText型に落としこむ際には気をつけなければならないことがあります。

以下の例は例外を投げてしまいます。

T.decodeUtf16BE "\xD8\x00"

出力

"*** Exception: Cannot decode input: Data.Text.Encoding.Fusion.streamUtf16BE: Invalid UTF-16BE stream

これは、Text型が内部的にはUTF-16を利用していることが原因です。 UnicodeではU+D800~U+DBFFとU+DC00~U+DFFFの範囲の16bitをペア(サロゲートペア)にして、16bitだけでは表現しきれない文字を表現しています。この範囲はUTF-16では表現しきれないので、decodeUtf16BEなどでは例外を投げます。外部から与えられたByteStringをデコードするときは、この挙動に注意して例外を拾うか、 decodeUtf16BEWith関数など、デコードのエラーが生じた時の挙動を明示する関数を使うといいでしょう。

lenientDecodeを使うと、デコードできないものはコードポイントU+FFFDに置き換えられます。

T.putStrLn $ T.decodeUtf16BEWith T.lenientDecode "\xD8\x00"

他にもUTF-8に限ればutf8-stringというパッケージがあるようですが、使ったことはありませんので分かりません(ヽ´ω`)

blaze-builder

チャンクのサイズが...

ByteStringの説明のところで、LazyなByteStringはStrictなByteStringのリストだと説明しました。だとすると、細かく作ったStrictなByteStringをたくさんつなげばLazyなByteStringを作る場合、結局細かいByteStringが連なったリストが出来上がりそうです。

下の例では、5バイトのByteStringを10個連結して、50バイトのLazyなByteStringを作っています。それをまたチャンクに分解して、それぞれの長さを調べると、確かに5バイトのチャンクが10個連なってます。

{-# LANGUAGE OverloadedStrings #-}

import Control.Monad
import Data.Monoid
import Blaze.ByteString.Builder           as B
import qualified Data.ByteString.Lazy     as BL
import qualified Data.ByteString.Char8    as BS


chunkSizes :: BL.ByteString -> [Int]
chunkSizes = map BS.length . BL.toChunks

concatByteString :: Int -> BS.ByteString -> BL.ByteString
concatByteString n = BL.fromChunks . replicate n

main = print . chunkSizes $ concatByteString 10 "xxxxx"

出力

[5,5,5,5,5,5,5,5,5,5]

5バイトずつのチャンクになって、これじゃァCharのリストだったStringと大して変わらんってことで、ひとつひとつのチャンクサイズを広げて、うまいことLazyなバイトストリングを作ってくれるのがblaze-builderです。

blaze-builderではBuilder型を介して文字列の結合をO(1)で行い、LazyなByteStringを吐き出すときに大きなチャンクを提供してくれます。

blaze-builderで文字列連結をするため、blazeConcatByteString関数を加えた例が以下のコードです。

{-# LANGUAGE OverloadedStrings #-}

import Control.Monad
import Data.Monoid
import Blaze.ByteString.Builder        as B
import qualified Data.ByteString.Lazy  as BL
import qualified Data.ByteString.Char8 as BS


byteN :: Int -> BS.ByteString
byteN n = BS.concat $ replicate n "x"

chunkSizes :: BL.ByteString -> [Int]
chunkSizes = map BS.length . BL.toChunks

concatByteString :: Int -> BS.ByteString -> BL.ByteString
concatByteString n = BL.fromChunks . replicate n

blazeConcatByteString :: Int -> BS.ByteString -> BL.ByteString
blazeConcatByteString n = B.toLazyByteString . mconcat . map B.fromByteString . replicate n

main =
    forM_ [1, 10, 100, 1000, 5000, 10000] $ \x -> do
        let bs = byteN x
        putStrLn $ show x ++ "-byte:"
        putStr "no-blaze "
        print . chunkSizes $ concatByteString      10 bs
        putStr "blaze    "
        print . chunkSizes $ blazeConcatByteString 10 bs
        putStr "\n"

出力

1-byte:
no-blaze [1,1,1,1,1,1,1,1,1,1]
blaze    [10]

10-byte:
no-blaze [10,10,10,10,10,10,10,10,10,10]
blaze    [100]

100-byte:
no-blaze [100,100,100,100,100,100,100,100,100,100]
blaze    [1000]

1000-byte:
no-blaze [1000,1000,1000,1000,1000,1000,1000,1000,1000,1000]
blaze    [4080,5920]

5000-byte:
no-blaze [5000,5000,5000,5000,5000,5000,5000,5000,5000,5000]
blaze    [4080,32752,13168]

10000-byte:
no-blaze [10000,10000,10000,10000,10000,10000,10000,10000,10000,10000]
blaze    [10000,10000,10000,10000,10000,10000,10000,10000,10000,10000]

blaze-builderを使ってつなげたほうが、チャンクのサイズが大きくなって、連結リストの長さが短くなってますね。こうしてチャンクのサイズを大きくすることで、キャッシュの利用を効率化したり、入出力に関わるシステムコールのオーバーヘッドを低減しているようです。

詳しいバッファ利用の戦略などについては、blaze-builderのhaddockや作者のブログ記事を参考にして下さい。

これを利用して、HTMLやSVGなどを高速に構築するblaze-htmlやblaze-svgといったパッケージがあります。 WebアプリケーションフレームワークのYesodなどでも使われていますね。

ベンチマーク

最後に、blaze-builderがどれくらい効くのが、ベンチマークをカジュアルに取ってみましょう！

ベンチマークにはcriterionを使ってます。Haskellのベンチマークが簡単に取れて、ブラウザで見られる出力も得られます。

blazeを使用しない元々の連結操作と、blazeを利用した連結との両方の方針で文字列片をつなげて、ファイルに出力しています。ついでに、Textに関してもByteStringと同様にベンチマークを取ってみました。

{-# LANGUAGE OverloadedStrings #-}

import Data.Monoid
import Blaze.ByteString.Builder           as B
import Blaze.ByteString.Builder.Char.Utf8 as BT
import qualified Data.ByteString.Lazy     as BL
import qualified Data.ByteString.Char8    as BS
import qualified Data.Text                as T
import qualified Data.Text.Lazy           as TL
import qualified Data.Text.Lazy.IO        as TL
import qualified Data.Text.Encoding       as T
import qualified Data.Text.Lazy.Builder   as T
import Criterion.Main
import System.IO (openTempFile)
import System.Directory (removeFile)


-- N-byteのByteStringを作る
byteN :: Int -> BS.ByteString
byteN n = BS.concat $ replicate n "x"

concatByteString :: Int -> BS.ByteString -> BL.ByteString
concatByteString n = BL.fromChunks . replicate n

blazeConcatByteString :: Int -> BS.ByteString -> BL.ByteString
blazeConcatByteString n = B.toLazyByteString . mconcat . map B.fromByteString . replicate n


-- N-文字のTextを作る
textN :: Int -> T.Text
textN n = T.concat $ replicate n "鬼"

concatText :: Int -> T.Text -> TL.Text
concatText n = TL.fromChunks . replicate n

blazeConcatText :: Int -> T.Text -> TL.Text
blazeConcatText n = T.toLazyText . mconcat . map T.fromText . replicate n


main :: IO ()
main = do
    -- 出力用の一時ファイル
    (tmpfile, tmph) <- openTempFile "/tmp" "bench_blaze"

    -- 1, 10, 50の大きさの単語を、10000個つなげる
    let samples = [1,5,10,50]
        wc      = 10000

    defaultMain $
        [ bgroup ("[bytestring] " ++ show x ++ "-byte words") [
              bench "no-blaze" $ BL.hPut tmph $ concatByteString      wc (byteN x)
            , bench "blaze"    $ BL.hPut tmph $ blazeConcatByteString wc (byteN x)
            ]
        | x <- samples ]
        ++
        [ bgroup ("[text] " ++ show x ++ "-letter words") [
              bench "no-blaze" $ TL.hPutStr tmph $ concatText      wc (textN x)
            , bench "blaze"    $ TL.hPutStr tmph $ blazeConcatText wc (textN x)
            ]
        | x <- samples ]

    removeFile tmpfile

上半分がByteStringで、下半分がTextです。

f:id:bicycle1885:20121224233049p:plain

ちっちゃい文字列片ほど、効果が大きいようですね。

ベンチマークの実行は以下の様な感じです。

% ghc -O2 bench_blaze.hs
% ./bench_blaze -o bench_blaze.html

出力は分かりやすさのため一部削ってます。

出力

warming up
estimating clock resolution...
mean is 1.488974 us (640001 iterations)
found 4828 outliers among 639999 samples (0.8%)
  3539 (0.6%) high severe
estimating cost of a clock call...
mean is 58.44394 ns (8 iterations)
found 1 outliers among 8 samples (12.5%)
  1 (12.5%) high severe

benchmarking [bytestring] 1-byte words/no-blaze
mean: 1.221033 ms, lb 1.214946 ms, ub 1.228679 ms, ci 0.950
std dev: 34.94729 us, lb 29.31052 us, ub 43.71415 us, ci 0.950

benchmarking [bytestring] 1-byte words/blaze
mean: 136.0696 us, lb 125.3679 us, ub 144.9379 us, ci 0.950
std dev: 49.58111 us, lb 39.05003 us, ub 60.22993 us, ci 0.950

benchmarking [bytestring] 5-byte words/no-blaze
mean: 1.321834 ms, lb 1.310183 ms, ub 1.337639 ms, ci 0.950
std dev: 68.97484 us, lb 54.09187 us, ub 95.68489 us, ci 0.950

benchmarking [bytestring] 5-byte words/blaze
mean: 743.6934 us, lb 709.6834 us, ub 775.6366 us, ci 0.950
std dev: 169.3327 us, lb 132.2535 us, ub 223.0841 us, ci 0.950

benchmarking [bytestring] 10-byte words/no-blaze
mean: 1.396986 ms, lb 1.387786 ms, ub 1.408378 ms, ci 0.950
std dev: 52.24961 us, lb 42.96288 us, ub 63.49543 us, ci 0.950

benchmarking [bytestring] 10-byte words/blaze
mean: 1.448711 ms, lb 1.383338 ms, ub 1.513390 ms, ci 0.950
std dev: 331.6004 us, lb 294.9020 us, ub 414.8301 us, ci 0.950

benchmarking [bytestring] 50-byte words/no-blaze
mean: 6.478227 ms, lb 5.424632 ms, ub 7.761446 ms, ci 0.950
std dev: 5.944880 ms, lb 5.028813 ms, ub 7.670362 ms, ci 0.950

benchmarking [bytestring] 50-byte words/blaze
mean: 7.454674 ms, lb 7.041077 ms, ub 8.173508 ms, ci 0.950
std dev: 2.719382 ms, lb 1.760055 ms, ub 5.017425 ms, ci 0.950

benchmarking [text] 1-letter words/no-blaze
mean: 3.269654 ms, lb 3.259670 ms, ub 3.282536 ms, ci 0.950
std dev: 58.06780 us, lb 46.70660 us, ub 75.66644 us, ci 0.950

benchmarking [text] 1-letter words/blaze
mean: 306.1035 us, lb 304.2076 us, ub 308.4018 us, ci 0.950
std dev: 10.65997 us, lb 9.045305 us, ub 12.66282 us, ci 0.950

benchmarking [text] 5-letter words/no-blaze
mean: 4.930164 ms, lb 4.542772 ms, ub 6.752361 ms, ci 0.950
std dev: 3.662762 ms, lb 177.2833 us, ub 8.698922 ms, ci 0.950

benchmarking [text] 5-letter words/blaze
mean: 1.524795 ms, lb 1.515778 ms, ub 1.535984 ms, ci 0.950
std dev: 51.46477 us, lb 43.33059 us, ub 64.69094 us, ci 0.950

benchmarking [text] 10-letter words/no-blaze
mean: 5.897089 ms, lb 5.879782 ms, ub 5.917257 ms, ci 0.950
std dev: 95.94202 us, lb 82.75690 us, ub 114.9103 us, ci 0.950

benchmarking [text] 10-letter words/blaze
mean: 3.982748 ms, lb 3.591212 ms, ub 4.620310 ms, ci 0.950
std dev: 2.510847 ms, lb 1.805720 ms, ub 3.784512 ms, ci 0.950

benchmarking [text] 50-letter words/no-blaze
mean: 21.43716 ms, lb 19.98548 ms, ub 27.01497 ms, ci 0.950
std dev: 13.20911 ms, lb 2.446692 ms, ub 31.00659 ms, ci 0.950

benchmarking [text] 50-letter words/blaze
mean: 21.46007 ms, lb 20.58942 ms, ub 22.80350 ms, ci 0.950
std dev: 5.445815 ms, lb 3.780456 ms, ub 8.059161 ms, ci 0.950

ところで、最新のbytestring-0.10.2.0のパッケージを見ると、Builder関係のパッケージがbytestringから提供されてますね。ドキュメントも、サンプルコードなどがあり充実してます。 textでは既にBuilderが入ってきてます。

まとめ

３種の文字列型と連結ライブラリのblaze-builderを見てきました。

Stringはマルチバイト文字も扱えてカジュアルに使えるけど長い文字列には向かない。
ByteStringは文字列というよりバイト列で、各種データをシリアライズするには便利。blaze-builderもあるし。
TextはStringのようにマルチバイト文字も扱えて、ByteStringのように効率のいい処理ができて嬉しい。

これらの実装もHaskellで書かれているので、ソースを読んでみればきっと面白い発見があると思いますよ(・ω<)。

何か間違いなどありましたら、Twitter宛@brain_appleかコメント欄などで教えて頂けるとありがたいです。

りんごがでている

何か役に立つことを書きます

JuliaTokyo #1 を開催しました

メインセッション

Juliaのこれまでとこれから - @sorami

Julia100本ノック - @chezou

Juliaのパッケージを作ろう! - @bicycle1885

ビジュアライズJulia - @nezuq

メカ女子将棋 - @kimrin

A Brief JuliaCon Report / Natural Language Processing with Julia - Pontus Stenetorpさん

LT

Juliaでちゃんとしたコードを書く - 関さん

Julia0.3でランダムフォレスト - @gepuro

Plotly Julia API - Yasunobu Igarashiさん

LightTableを使いましょう - @QuantixResearch

感想

HaskellのFusionがあれば速度と抽象化を両立できる

ベクトル化された記法

分散の計算

fusionを使って宣言的に書き、効率的なコードを吐く

fusionがどうやって動いているのか

分散の計算、再び

パフォーマンス比較

新しいfusion

Generalized Stream Fusion

Data Flow Fusion

Pythonのコマンドライン引数処理の決定版 docopt (あとJuliaに移植したよ)

Julia版のdocopt - DocOpt.jl

Haskellでデータ処理がしたい！

Juliaのデストラクター

matplotlib入門

Haskell Tips (文字列編)

文字列

String

ByteString

Text

blaze-builder

チャンクのサイズが...

ベンチマーク

まとめ