プログラミングの実験場

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

HaskellでASTからのコード生成（１）　Writer+Stateモナドによる構造化

前のポストでGADTを使って型付きの抽象構文木（AST）を表現する方法について書いたが、ここではそのASTから他の言語のコード生成する方法について調査・検討した結果を記す。@keigoiさんの記事（http://d.hatena.ne.jp/keigoi/20111206/haskell_tagless_dsl）に触発され、一部借用したものである。なお、今回のコード生成の方法論はASTの型付き・型なしの区別にかかわらず使える。

MatlabのDSL in Haskell

今回例として扱うのは、Matlab（行列計算言語・統合環境）の行列計算を表現するようなDSL。最小構成として、２次元行列の加減乗算と要素の取り出しを表現したい。（Matlabの文法の簡単な説明はこれなどを参照。 http://www.math.meiji.ac.jp/~mk/labo/text/matlab/node4.html）これだけの小さなものでも、

  • 変数の管理
  • let-sharing（別の記事で記述予定）

といった、コード生成のエッセンスが含まれている。

まず、構文木を表すデータ型Exp rはこのような感じになる。

import Data.Text (Text)
import qualified Data.Text as T

type Ident = Text

-- 2D matrix
data Mat2D r = Mat Ident | MatConst [[r]]

class Scalar a
instance Scalar Int
instance Scalar Integer
instance Scalar Double

data Exp r where
	Const :: r -> Exp r
	Add :: Exp r -> Exp r -> Exp r
	Sub :: Exp r -> Exp r -> Exp r
	Mul :: Exp r -> Exp r -> Exp r
	Abs :: (Scalar r) => Exp r -> Exp r
	Signum :: (Scalar r) => Exp r -> Exp r
	Index2D :: (Num r) => Exp Int -> Exp Int -> Exp (Mat2D r) -> Exp r

この構文木を使って、たとえば、reify :: Exp r -> Textという関数を定義して、

> reify $ Const (MatConst [[1,2,3],[4,5,6]])
[[1 2 3];[4 5 6]]
> reify $ Index2D 1 2 (MatConst [[1,2,3],[4,5,6]])
m = [[1 2 3];[4 5 6]];
m(1,2)

といった具合のコード生成をしたい。「コード生成できるもの」を表すReifiableクラスというものを作って、reify :: a -> TextをReifiableクラスのメソッドとする。ここでReifiableのインスタンスとなるのは、Exp r, Mat2D r, Int, Integer, Doubleである。例として以下にExp rとMat2D rのインスタンス定義を示す。ここまでの全ソースコードと実行例はここにある。

instance Reifiable (Exp r) where
	reify (Const exp) = reify exp
	reify (Add e1 e2) = T.concat [reify e1,"+",reify e2]
	reify (Index2D i j mat) = T.concat [reify mat,"(",reify i,",",reify j,")"]

instance (Show r) => Reifiable (Mat2D r) where
	reify (Mat ident) = ident
	reify (MatConst xss) = T.pack $ "[" ++ intercalate ";" (map (intercalate " " . map show) xss) ++ "]"

ちなみに、ここでは行列は二次元のみであるが、行列を２次元以外にも対応させて型安全にするには、ここ（http://www.haskell.org/haskellwiki/GHC/Type_families）にあるように、data family Mat dim rなどとして、dimにD1, D2, D3といった型のタグを付けてやれば良いと思われる。（ただしそのようなことをやると、型クラスをつけるときにだいぶ複雑になる。）

ASTは単なるHaskellの値であるので、mapやfoldなどを使ってASTを構築することもできる。

*Main> reify $ foldl1 (+) $ replicate 5 $ Const (MatConst [[1,2,3],[4,5,6]])
"[1 2 3;4 5 6]+[1 2 3;4 5 6]+[1 2 3;4 5 6]+[1 2 3;4 5 6]+[1 2 3;4 5 6]"

なかなか良い感じなのだが、実はこの実装では問題があり、Index2Dによる要素の取り出しがきちんと出力できない。
Index2D 1 2 (MatConst [[1,2,3],[4,5,6]])をreifyすると、[[1,2,3],[4,5,6]](1,2)となるのだが、Matlabのいけてない言語仕様により、このようには書けず、定数行列を一回変数に代入する必要がある！

モナドを利用した構造化

  • そういった、コード出力時に単純変換以上の何らかの操作をしたい場合には、状態管理を出力のロジックに組み入れる必要がある。→ Stateモナドを使う。
  • また、コードの出力も、単にTextを返すのではなく、ASTをたどっていく過程の任意の段階で出力をしたい。→ Writerモナドを使う。

State、Writerの両方を使いたいので、モナド変換子で２つのモナドを積み上げるのが基本のようだが、今回はRWSTというモナド変換子（Hackage、Shelarcyさんによる解説）を使う。RWSモナドはReader, Writer, Stateのすべてが組み込まれており、しかもMonadState(http://hackage.haskell.org/packages/archive/mtl/2.0.1.0/doc/html/Control-Monad-State-Class.html)などの型クラスのインスタンスになっているので、get/put, tellなどがliftせずに使えるようになっている便利なモナドである。IOも、後ほどのlet-sharingの件で必要となるので入れている。

type W a = RWST () Text CompilerState IO a

最初この巨大な型を見ると面食らうが、その意味するところは、型パラメータの左から順に、読み込む（Reader）ものは無く（()）、Textを書き出すことができ、CompilerStateという状態を書き換えることができ、IOもliftIOを介して使える、というものである。書きだすもの（Text）と状態（CompilerState）は型が決まっており、get,put,tellなどの関数を使って、モナドの裏側で暗黙に操作される（これらの関数はみなW a型の値を返す。）。 型コンストラクタの最後に付くa型はdo記法の中で自由に替えられ、モナドのユーザーが明示的に扱うものであるというのは、IOなど他のモナドと同じである。

最初の、reifyがTextを返すバージョンの抜粋

class Reifiable a where
	reify :: a -> Text

instance Reifiable (Exp r) where
	reify (Const exp) = reify exp
	reify (Add e1 e2) = T.concat [reify e1,"+",reify e2]
	reify (Index2D i j mat) = T.concat [reify mat,"(",reify i,",",reify j,")"]

gen :: Exp r -> IO ()
gen = T.putStrLn . reify

reifyがモナドを返すように書き換えたバージョンの抜粋。

type W a = RWST () Text CompilerState IO a
class Reifiable a where
	reify :: a -> W Text

instance Reifiable (Exp r) where
	reify (Const exp) = reify exp
	reify (Add e1 e2) = do
		c1 <- reify e1
		c2 <- reify e2
		return $ T.concat [c1,"+",c2]
	reify (Index2D i j mat) = do
		ci <- reify i
		cj <- reify j
		m <- reify mat
		return $ T.concat [m,"(",ci,",",cj,")"]

data CompilerState = CompilerState { variableCount :: Int }

gen :: Exp r -> IO ()
gen e = do
	(txt,_,_) <- runRWST (reify e) () (CompilerState 0)
	T.putStrLn txt

まだ機能は古いバージョンと全く同じ。runRWSTを使ってモナドを走らせると、IO (Text, CompilerState,Text)を返すので、その結果を取り出せる。

先に述べたように、定数行列から要素を取り出すとき、いったん行列を変数に代入して、その変数から要素取り出しをしたい。
どのように設計すればよいかというと、

  • 以前と同様、reifyが返すのはその値を参照するのに必要なコード。つまり代入した変数の名前。
  • それに加えて、MatConstに対するreifyの定義の中で、tellを使って、値を準備するコード（つまり新しい変数への代入のコード）をWriterモナドに「出力」する。
instance (Show r) => Reifiable (Mat2D r) where
	reify (Mat ident) = return ident
	reify (MatConst xss) = do
		let code =  T.pack $ "[" ++ intercalate ";" (map (intercalate " " . map show) xss) ++ "]"
		var <- newName
		tell $ T.concat [var," = ",code,";\n"]
		return var

ここで、「出力」と鍵括弧にくくったのは、putStrLnのようにIOに出力するのではなく、あくまでWriterモナドをつかっているので、最後にモナドを走らせて「出力」を取り出す必要があるからである。newNameというのは、Wモナドが持つ状態（WモナドはMonadStateのインスタンスである）を利用して一意な名前を作る関数で、次のように定義できる。

data CompilerState = CompilerState { variableCount :: Int }

newName :: W Text
newName = do
	CompilerState count <- get
	put $ CompilerState (count + 1)
	return $ T.pack $ "v" ++ show (count+1)

最後に、Writerモナドに出力した内容と、reify自体の返り値をまとめてコード出力の結果とする。

gen :: Exp r -> IO ()
gen e = do
	(txt,_,assigned) <- runRWST (reify e) () def
	T.putStrLn $ T.append assigned txt

実行結果は以下のようになる。

*Main> let m1 = Const (MatConst [[1,2,3],[4,5,6]]) :: (Exp (Mat2D Int))
*Main> let m2 = Const (MatConst [[10,5,3],[1,4,2]]) :: (Exp (Mat2D Int))
*Main> gen $ Index2D (Index2D 1 2 m1) 1 m2
v1 = [1 2 3;4 5 6];
v2 = [10 5 3;1 4 2];
v2(v1(1,2),1)

良い感じに変数への代入ができている！　コードをここに挙げておく。

変数の重複の排除

ここで、以下の様な式を考えてみる。

m1 :: Exp (Mat2D Int)
m1 = Const (MatConst [[1,2,3],[4,5,6]])

test4 :: Exp (Mat2D Int)
test4 = m1 + m1

これをreifyすると、

*Main> gen test4
v1 = [1 2 3;4 5 6];
v2 = [1 2 3;4 5 6];
v1+v2

となり、同じ内容であるはずのm1がv1,v2と複数回生成されている。これくらいの小さなデータならいいが、もしこのm1が重い計算の結果得られる値の場合このような重複は避けたい。なぜこういうことが起こるかというと、test4が簡約されるときに、m1 + m1という式の両側のm1が同じがどうかというのを構わずにHaskellのコンパイラが展開してしまうからである。この論文の3.1節に説明されている。これの重複を避けることには、Let-sharingという名前が付いている。一言で言えば、System.Mem.StableNameモジュールのmakeStableNameとhashStableNameという関数を使って、項の同一性判定をしてやれば良く、さほど難しくはなかった。別の記事で実装について記す予定。

参考になった記事

文献を調査中なので、参考になるものがあったら教えてほしいです。

GADTによるHaskellの型付きDSLの構築

GADTについて、この解説が分かりやすい（英語）。 http://en.wikibooks.org/wiki/Haskell/GADT

背景

Haskell（というか関数型言語一般）の素晴らしい機能の一つ、代数的データ型を使うと、抽象構文木(AST)をHaskellの枠組みの中で表現できる。
例えば、整数の足し算・引き算ならば、

data Exp = Const Int | Add Exp Exp | Subtract Exp Exp

さらにこのASTを数値一般に使いたいので型パラメータrをつけて多相にしてみる。

data Exp r = Const r | Add (Exp r) (Exp r) | Subtract (Exp r) (Exp r)

これを使えば、例えば(3 + 4) - 2という数式は、Subtract (Add (Const 3) (Const 4)) (Const 2) :: Num r => Exp rというHaskellの（多相の）値として表現できる。Exp r型は、rという型パラメータによって型安全性が保たれていて、Add (Const (3::Int)) (Const (4::Double))というのはコンパイルが通らない。

いったんHaskellの値として表現出来れば、あとはパターンマッチを用いるなどしてHaskellの枠組みで型安全に、自由に操作できる。たとえばこのASTを元に、コードの評価をしたり、他の言語の表現として出力したり、といった様々なことができる。つまり、Haskellの中に小さな言語を構築し、Haskellをホスト言語・処理系として使う、一種のDSLといえる。

たとえば評価ならば、

eval :: (Num r) => Exp r -> r
eval (Const v) = v
eval (Add e1 e2) = (eval e1) + (eval e2)
eval (Subtract e1 e2) = (eval e1) - (eval e2)

といった具合。

コード生成だったら、

reify :: (Show r) => Exp r -> String
reify (Const v) = show v
reify (Add e1 e2) = concat ["(",reify e1,"+",reify e2,")"]
reify (Subtract e1 e2) = concat ["(",reify e1,"-",reify e2,")"]

もちろん、掛け算・割り算にも、データコンストラクタとそれに対応する実装をeval,reifyに追加することで簡単に対応できる。

ここまでのコードと実行例は： https://gist.github.com/nebuta/6096315

ASTを型付きにする。

ここで、先のDSLを拡張して、式の評価結果を比較するEq演算子を使えるようにしてみたい。

data Exp r = Const r | Add (Exp r) (Exp r) | Subtract (Exp r) (Exp r)
			| Eq (Exp r) (Exp r) -- ??

ここで問題が発生する。上のような通常のデータ型の定義だと、すべてのデータコンストラクタは引数を適用するとExp r型になり、このEqデータコンストラクタが返す型は、たとえば、(Eq (Const (1::Int)) (Const 1)) :: Exp IntとようにEqのオペランドと同じ型になってしまう。欲しいのはExp IntではなくてExp Boolである。

そこで使えるのがGADTである。GADT(Generalized algebraic datatype)とは、多相のデータ型に対し、データコンストラクタの型を特化することができる仕組み（というのが私の理解）。
まずは元のExp rの定義をGADTのスタイルに書き換える。

{-# LANGUAGE GADTs #-} -- GADTsのプラグマが必要。
data Exp r where
	Const :: r -> Exp r
	Add :: Exp r -> Exp r -> Exp r
	Subtract :: Exp r -> Exp r -> Exp r
	Eq :: Exp r -> Exp r -> Exp r

この時点では、前の例と意味は全く同じ。データコンストラクタが関数の型シグネチャのスタイルに書換えられているのが分かる。型シグネチャの右端に現れるのが、データコンストラクタに引数を完全に適用した際に得られる型。つまりどれもExp r型である。

今やりたいのは、Eq (Exp r) (Exp r)の型をExp rではなくて、より特化したExp Boolにしたいということ。GADTでは、このデータコンストラクタが返す値の型を、Exp rを特化した任意の型に置き換えることができる。

	Eq :: Exp r -> Exp r -> Exp Bool

に書き換えれば良い。

data Exp r where
	Const :: r -> Exp r
	Add :: Exp r -> Exp r -> Exp r
	Subtract :: Exp r -> Exp r -> Exp r
	Eq :: Exp r -> Exp r -> Exp Bool

さらに、２つのExp Boolを引数にとるような、AndとOrも定義してみる。

	And :: Exp Bool -> Exp Bool -> Exp Bool
	Or :: Exp Bool -> Exp Bool -> Exp Bool

GADTを使わないバージョンではrは数値型というのを仮定していたが、この新しいバージョンではrはBoolもありうる。だが、AddやSubtractの引数にはExp Boolは取りたくない。そこで、以下のようにAdd, Subtractの型パラメータrに制限をつけてやる。

data Exp r where
	Const :: (Show r) => r -> Exp r
	Add :: (Num r) => Exp r -> Exp r -> Exp r
	Subtract :: (Num r) => Exp r -> Exp r -> Exp r
	Eq :: (Eq r) => Exp r -> Exp r -> Exp Bool
	And :: Exp Bool -> Exp Bool -> Exp Bool
	Or :: Exp Bool -> Exp Bool -> Exp Bool

こういうことができるのもGADTの強み。
そして、evalの実装は以下のように変わる。

eval :: Exp r -> r
eval (Eq e1 e2) = eval e1 == eval e2
eval (And e1 e2) = eval e1 && eval e2
eval (Or e1 e2) = eval e1 || eval e2
{- 他のパターンは前と同じ -}

evalの定義に(Num r)の制限がなくなったことに注意。これは、型クラスの制限がデータコンストラクタについたために、関数定義には制限をつける必要がなくなったからである。

これによって、And (Eq (Const 3) (Const (2::Int))) (Const True) :: Exp Bool のような、内部に複数の型が混ざったような構文木が型安全に表現できる。内側にあるConst (2::Int)の型はExp Intだが、最終的にAndを適用した後はExp Boolになる。Exp rという多相型の枠組みの中で、型を柔軟かつ強力に付けられる利点がある。

GADTに書き換えたコードと実行例は：https://gist.github.com/nebuta/6096345

コード生成

この型付きASTから手続き型言語のコードを生成する方法についても調査・検討したので別の記事で書く予定。この記事にあるように、WriterモナドとStateモナドを組み合わせることで、一意な変数名の生成などの複雑な処理が書ける。
http://d.hatena.ne.jp/keigoi/20111206/haskell_tagless_dsl

この手法を使うことで、画像処理、JavaScript、いろいろな型安全DSLをHaskell内で構築できると思われ、いろいろなアイデアが浮かんでくる。

また、この論文（Haskell Symposium 2010および2009）はlet-sharingとlambda-sharingという、HaskellのDSLからのコード生成についての２つの問題についての解決法がわかりやすく書いてある。
http://www.eecs.harvard.edu/~mainland/publications/mainland10nikola.pdf
http://www.cs.uu.nl/wiki/pub/Afp/CourseLiterature/Gill-09-TypeSafeReification.pdf

D3.jsを使ってブラウザ上で宣言的なデータ可視化

D3.jsというJavaScriptでのDOM操作・可視化ライブラリ http://d3js.org/ が面白そう。D3.jsはHTMLのDOM要素とSVGに対して宣言的にプログラミングすることで可視化が見通しよく書けるライブラリ。D3.jsで生成されるのは図形を含めてすべてDOM要素であるので、CSSや他のJSコードと自在に組み合わせられるという利点がある。便利関数がたくさんあって多機能な感じ。jQueryのような感覚でDOM要素の集合を選択し、それに対してデータ（配列、ハッシュの配列、配列の配列）を割り当てて、DOM要素の属性を変更できる。

資料

最小構成

<html>
<head>
	<title>D3.js test</title>
</head>
<body>
	<div id='div1'></div>
	<div id='div2'></div>

	<!-- 必要なJSはこのファイル１つだけ。JQueryなども不要。 -->
	<script src="http://d3js.org/d3.v3.min.js" charset="utf-8"></script>
	<!-- 自分のコード -->
	<script src="d3test.js" charset="utf-8"></script>

</body>
</html>

DOM操作のライブラリなので、JSファイル１つのみで、専用のCSSなどは無い。

基本的なDOM操作

jQueryのように、DOM要素を選択し、メソッドチェーンで子要素を追加したり、属性や内容を書き換えられる。

function test0(){
	d3.select('body').append('div').attr('id','graph').text('Hey');
}

SVG要素を操作する例。

function test1(){
	//SVG要素を作成。
	var svg1 = d3.select("#div1").append("svg")
		.attr("width",300)
		.attr("height",300)
		.style("background","#eef");

	//子要素を追加。
	svg1.append('circle')
		.attr("cx", 100)
	    .attr("cy", function() {return 10+Math.random()*200;}) //値の指定に、関数で値を返すこともできる。
	    .attr("r", 20)
	    .style('fill',function(){return 'pink';})
		.on("mouseover", function(){d3.select(this).style("fill", "gray");})    //イベントハンドラの指定。（この動作ならCSSで書いたほうが良さそうだが。）
	    .on("mouseout", function(){d3.select(this).style("fill", "pink");});
}

２次元配列からデータをとってきて、テーブルのセルを生成する例

//参考：http://christopheviau.com/d3_tutorial/
function test2(){
    var array2d = [[1,3,2],[4,9,6],[1,2,4]];

    // selectAllでデータのプレースホルダを作り、それにdataでデータを追加する。そしてenterでデータに対応する集合が実際に出来る。
    // 解説はこれが分かりやすい。http://bost.ocks.org/mike/join/

    d3.select("#div2 table")
		.style("background","#fee")

	    .selectAll("tr")　　　//selectAll,data.enterという一組でデータを追加。
	    .data(array2d)
	    .enter().append("tr")	//データarray2dの各要素（array2d[0],array2d[1],array2d[2]）に対応して行を作成。
	    
	    .selectAll("td")	//各行に対して、列の作成
	    .data(function(d){	
	    	console.log(d);
	    	return d;
	    })
	    .enter().append("td")

	    .style("background",function(d,i,j){ //dは２次元配列の中身、(i,j)は２次元配列のインデックス。
	    	var d_max = 9;
	    	var d_min = 0;
	    	console.log(i,j);
	    	//HSLカラーでヒートマップを表現。
	    	//色の指定については、https://github.com/mbostock/d3/wiki/Colors
	       	return d3.hsl(240*(1-d/d_max),1,0.5);
	    })
	    .attr("width",50)	//各tdのサイズを指定。
	    .attr("height",50)

	    .text(function(d,i,j){ return "("+i+","+j+"): "+d; });
}

一番最後の例など、メソッドチェーンが長大。このような入れ子の構造でメソッドの返り値が何の「型」なのか、そして選択されたDOMの集合がどう変わるのか、というのがまだちょっと理解できていない。
こういうところこそ、Haskellのような型が便利になる場面なのだと思う。Haskellのライブラリは、型があるので理解が早い。と、無理やりHaskellにつなげた所で今日のところはここで終わり。

MIT Halideのインストール方法

C++（とOCaml）で書かれた高速画像処理ライブラリである MIT Halide の自分の環境（Mac OSX 10.8）でのインストール方法の覚書。

ここからコンパイル済みのものがダウンロードできる。解説のスライドや論文へのリンクもここに載っている。
MIT Halide
http://halide-lang.org/

アーカイブを展開するとlibHalide.aとHalide.hという２つのファイルができるので、然るべき場所（/opt/local/includeと/opt/local/libなど）にコピーしておく。

プログラムの例はHalideのgitレポジトリを持ってくれば、testフォルダとappsフォルダに入っている。コンパイルは以下のようにg++を使ってやればよい。僕の環境にはg++が、g++-4.2, g++-mp-4.6, g++-mp-4.7 といろいろ入っているのだが、バージョン>=4.6が必要と思われる。

$ git clone https://github.com/halide/Halide
$ cd Halide/test
$ g++-mp-4.6 -std=c++0x vectorize.cpp -L halide -ldl -lpthread -I/opt/local/include /opt/local/lib/libHalide.a
$ ./a.out
HL_TARGET not set - inferring from host architecture...got x86_64
1.37134e+12 1.37134e+12 1.37134e+12
Vectorized vs scalar (float x 4): 0.519ms 1.76ms. Speedup = 3.400
1.37134e+12 1.37134e+12 1.37134e+12
Vectorized vs scalar (float x 8): 0.819ms 3.41ms. Speedup = 4.167
1.37134e+12 1.37134e+12 1.37134e+12
Vectorized vs scalar (double x 2): 0.687ms 0.917ms. Speedup = 1.335
1.37134e+12 1.37134e+12 1.37134e+12
Vectorized vs scalar (uint8_t x 16): 0.515ms 5.69ms. Speedup = 11.045
1.37134e+12 1.37134e+12 1.37134e+12
Vectorized vs scalar (int8_t x 16): 0.587ms 8.41ms. Speedup = 14.334
1.37134e+12 1.37134e+12 1.37134e+12
Vectorized vs scalar (uint16_t x 8): 0.551ms 3.02ms. Speedup = 5.477
1.37134e+12 1.37134e+12 1.37134e+12
Vectorized vs scalar (int16_t x 8): 0.523ms 3.41ms. Speedup = 6.521
1.37134e+12 1.37134e+12 1.37134e+12
Vectorized vs scalar (uint32_t x 4): 0.516ms 1.5ms. Speedup = 2.914
1.37134e+12 1.37134e+12 1.37134e+12
Vectorized vs scalar (int32_t x 4): 0.553ms 1.43ms. Speedup = 2.584
Success!

なお、まだAPIの仕様変更があるためか、公式サイトから持ってきたバイナリでは、Githubの最新のテストコードの中には、関数が見つからないといってコンパイルできないものが多かった。

追記：https://github.com/halide/Halide/blob/master/INSTALLにしたがってビルドできた。LLVMの3.2以上が必要（sudo port install llvm-3.2で入る）。それと、僕の環境では、llvm-configをln -s /opt/local/bin/llvm-config-mp-3.2 /opt/local/bin/llvm-configで作ってやる必要があった。makeでbin/libHalide.a, bin/libHalide.so, と include/Halide.hが生成され、しかるべきところにコピーすれば使える。ちゃんと最新のtestフォルダ内のコードもコンパイルが通った。

関数型言語による画像処理環境の追求　その２　Haskellで画像データを扱うためのアプローチ

Haskellで画像データを扱う際に、いろいろな手法がありうる。

Pure Haskell

Hackageを見てみると、（denseな）行列を扱うpure Haskellの（＝外部ライブラリに依存しない）ライブラリで有名なのは、これだろうか。並列化が自動でできるというのが売りのようだ。Simon Peyton JonesもYouTubeに載っている講演で紹介していたと思う。

shelarcyさんによる解説 http://itpro.nikkeibp.co.jp/article/COLUMN/20120605/400424/

ただ、速度は遅く、ちょっと実用には辛いレベルかと思われる。 http://repa.ouroborus.net/wiki/Examples/Fft2dHighpass

やはり、Pure Haskellで大きな数値データを高速に扱うのは未だ難しい（というか、あまり真剣に検討されていない？）という印象だ。

FFIで呼び出し

別にpure Haskellに拘る必要はなく、Cで書かれた外部ライブラリを呼び出すバインディングでも良い。と思って今一度Hackageを探してみたのだが、まともなのは殆ど見つからなかった。OpenCVへのバインディングが２つあるが、以下のものはわりと型もちゃんとついて良い感じに見える。

http://new-hackage.haskell.org/package/CV

それと、画像に特化せず行列のライブラリであれば、以下のようなものがある（GSL, LAPACK, BLASのバインディング）。

http://new-hackage.haskell.org/package/hmatrix 

ただこれらも、ベストではないように思える。というのも、immutableなデータが多く生成されてしまう可能性がある。例えばフィルターを掛けるたびに一時データの画像が生成・コピーされ、メモリとCPUの効率が悪い。FFIを繰り返し呼び出すとそれがパフォーマンス上のネックになる可能性もある。ちなみに、この問題を解決するためには、中間データの生成を無くす方法である、GHCの書換え規則や融合変換などが使えそう。融合変換ってリストだけでなく任意のデータ構造に使えるのだろうか。

（融合変換についての解説：http://itpro.nikkeibp.co.jp/article/COLUMN/20100706/349960/ http://big.freett.com/shelarcy/rewriterules.pdf）

HaskellからCのコードを生成

Haskellのプログラムにembedded DSL (EDSL)を埋め込み、EDSLとして記述された画像処理のアルゴリズムからC（あるいはCUDAなど）のソースを出力する方法。２種類が考えられる。

  1. Haskellプログラムの実行時にCソース出力し、それをCコンパイラでコンパイルする。
  2. Template Haskellを使ってコンパイル時にCソースを出力し、それをHaskellからリンクする。

Haskell側でのアルゴリズムの最適化なども可能で、実現できるならこれが一番優れているように思える。ただし、Hackageを見た＋軽くウェブ検索した限り、このようなアプローチの実装は殆ど無いようだ。

数少ない例としては、Nikolaという、Haskell内の embedded DSLからCUDAのコードを生成するライブラリがある。上でいう2の方法に当たる。論文が出たのが2010年と結構最近のものなので、たぶんこのような例は多くないことが推測される。

http://www.eecs.harvard.edu/~mainland/projects/nikola/

http://www.eecs.harvard.edu/~mainland/publications/mainland10nikola.pdf 論文

また、画像処理ではないが、Paraisoというライブラリは、Haskellでアルゴリズムを書くと、微分方程式を解くCUDAあるいはC+OpenMPのコードを生成する。これは上でいう1の方法にあたる。Hackage初出が2011年。

http://new-hackage.haskell.org/package/Paraiso

 

まとめ＋今度の展望

HaskellのライブラリのすべてがHackageに載っているわけではないので、どこか見落としている可能性はあるが、現状で普通のユーザーに使える本格的なHaskellの画像処理ライブラリはほとんどない、というのが現時点での調査の結論。

これから何らかの画像処理ライブラリを自分で作るとしたら、Haskell内のEDSLからCコードを生成するアプローチが一番有望（大変そうではあるが）のように思えるので、このアプローチについてこれからヒマを見て調べてみようと思う。

Haskellで型付きのDSL（型付きじゃないとHaskellに拘る意味が無くなってしまうので）を作る方法については、以下の資料が参考になりそう。

http://d.hatena.ne.jp/keigoi/20111206/haskell_tagless_dsl

http://okmij.org/ftp/tagless-final/index.html

 

関数型言語による画像処理環境の追求　その１　現状と背景

Haskellグループができたのを機にはてなブログを作ったので、せっかくなのでHaskellについてなにか書こうと、僕の専門分野で使う、Haskellをつかった生物画像処理の環境について今一度調べてみることにした。

生物研究での画像処理というのは主に、画像に対して処理をして特徴抽出したり、物体検出（対象は細胞や蛍光分子）してその位置・数・大きさ・形などを数えたり、分類したりする、ということである。 解析する画像データは通常denseな2〜5次元（平面 × Z軸 × チャンネル（異なる色の蛍光分子） × 時間）の行列であり、大きさは面あたり数百×数百、Z軸は数１０、チャンネルは1-5、時間軸は1-1000個程度のものを使うことが多い。

以下のリンクが参考になる。

Image analysis for biology http://www.rpgroup.caltech.edu/courses/PBoC%20ASTAR/files_2011/Matlab/Image%20Analysis%20with%20Matlab.pdf

この2010年の論文（PLoS Computational Biology）なんかもイントロとして良さそう。

Pattern Recognition Software and Techniques for Biological Image Analysis: http://www.ploscompbiol.org/article/info%3Adoi%2F10.1371%2Fjournal.pcbi.1000974 

 

これまで僕は画像処理について、MATLABで行き当たりばったりでやっていたのだが、MATLABにはメリットとデメリットがある。 まずはこの記事で、その現状について少し書いてみる（なので、これはまだHaskellの記事ではない）。

MATLABのメリット

  • 組み込み関数はそこそこ速い。
  • 使用例がそこそこある。ウェブで見つかる大学の講義資料などはたいてい（学習しやすいという理由で）MATLABを使っている。また公式サイトのドキュメントも詳しく丁寧。
  • オールインワンで、GUIでの出力も簡単。デバッガが使え、ステップ実行や、変数内容の表示や変更ができる。
  • 生物系の画像処理ではシェアがかなり高い印象。

デメリット

  • 言語仕様やライブラリが主流のモダン言語と比べて圧倒的に劣っている。仕様がチグハグで一貫性がなく、抽象化能力も低い。バッチ処理で必要になる、文字列・ファイル処理も不便。プログラムの再利用性が低い。
  • 変数に型がない。インタプリタなので、実行前の静的チェックができない。
  • 組み込み関数を使わないで行列を手動で操作したりすると実行速度が遅い。
  • ソフトが商用で、非常に高い。アカデミックライセンスもライセンスチェックが厳しい。
  • 組み込みのエディタがしょぼい。

このような考察の結果、もうちょっとクリーンで快適かつ型安全な画像処理プログラミングをしてみたいと思った次第。抽象化能力の高いモダンな言語、できればHaskellで出来ればいいのだけど、ウェブ上でほとんど実用例は見つからない。

ちなみに、MATLAB以外でのオプションは以下の様な感じだろうか。

  • C++ + OpenCV（あるいは他のライブラリ）
  • ImageJ (or Fiji) + プラグイン

ImageJはJavaで書かれたオープンソースで、Java（あるいはJVM上の言語、Scalaなど）で書いたプラグインで拡張可能。実際僕もImageJのプラグインをScalaで書いて、それで画像解析をした。Fijiはそれにプラグインを追加したもので、JRubyやclojureなどでスクリプトを書ける。ImageJはAPIの筋がいいし、Scalaが使えるので悪くない。

次のポストで、Haskellでの画像処理、行列の扱いなどについて現時点でのライブラリの状況についてまとめる。

もちろん、道具選び・道具作りに溺れて実際の解析にとりかかるまでの準備に時間をかけすぎては本末転倒なのだが、もっと計算の意味を高い抽象化で表現して、実装の詳細と分離するようなことができればプログラムの見通しも良くなるし、生産性が上がることは間違いない。

たとえば、MITがHalideというC++用の画像処理ライブラリを最近公開したが、これは純粋関数的な「アルゴリズム」と、実装の詳細の「スケジュール」を分けることで最適化を容易にしている。

http://people.csail.mit.edu/jrk/halide12/

この方向性は関数型言語の考えとも共通するし、並列化などを考えても非常に有望だと思う。