Haskellの配列周辺についてのメモ（続々）

Repaライブラリについて、簡単な備忘録をつけておく。

なにはともあれこのチュートリアルを読むわけだが、自分のために以下で簡単にまとめておく。

• repa (Data.Array.Repa) 基本となるライブラリ
• repa-io (Data.Array.Repa.IO) 画像を含む配列データのファイル入出力支援
• repa-algprithms (Data.Array.Repa.Algorthms) FFTとかDFTとかのライブラリ

Repa配列の型とかについて

インデックス

Repa配列のインデックスはやや特殊で、IArrayやMArrayのようなIx型クラスではなくShape型クラスを使う。Shape型クラスのインスタンスは特殊で、Zというゼロ次元配列を表す型か、これに中置型構築子:.を使ってInt型を必要な次元数連ねて構築されるZ :. Int :. Int :. Intのような型になる。なお:.は左結合なのでこれは(((Z:.Int):.Int):.Int)としても同じ（ドキュメントではこの書き方がよく見られる）。これは表記が面倒なので、Zに付加されたIntの個数（つまり配列の次元数）に応じて、DIM1 = Z:.Int, DIM2=Z:.Int:.Int, ...といった型シノニムが DIM5まである（とはいえ型シノニムなので型注釈でしか使えないし、どうせ同じものをデータ構築子としてパターンマッチで頻繁に使うことになるのでそこまで嬉しくない）。Zと(:.)は値構築子としても使うので、Shape型クラスのインスタンスの値はたとえば(Z:.99)::DIM1みたいになる。

配列の種類

配列の要素に何が入っているかによって、Repaの配列は何種類かに分けられる。たとえば：

• U (Unboxed) : 即値
• D (Delayed) : インデックスから値への関数（要はサンク値）

Dは重要で、たとえばU型のRepa配列に対してmapとかなにかの計算を実行すると、返ってくる値はD型になり、実際の計算が遅延される。D型のRepa配列は格納されている計算を実行することでU型にすることができる。この計算をライブラリの名前通り並列的 pararellに実行するのがComputeP函数。代わりに逐次的 sequentialに実行するのがComputeS函数。ComputeP函数は注意が必要で、並列性の関係で返ってくる値がモナドに包まれて返ってくる。もちろん純粋な函数の計算に副作用があるわけではなく、STモナドの中で使ってrunSTで普通に取り出せばなにも問題ない：

addArray f g = runST $ do {computeP $ zipWith (+) f g;}
h = addArray f g

こんな感じ。しかし、実はどんなモナドでも構わないので、リストモナドで：

[h] = computeP $ zipWith (+) f g

とかしても、STモナドの場合とまったく同じように、計算して即値になった配列がhに束縛されて返ってくる。まあでもSTモナドが一番わかり易いと思うけど。IOモナドでもいいのだが、そうするとモナド外に取り出せなくなるのでIOモナドの中で配列を操作しているとき以外にはオススメできない（基本は別のモナドを使ってからそのモナドをひん剥いて、配列を操作する函数全体はComputeSの場合と同じようにあくまで通常の純粋な函数にしておくべきだと思う）。とにかくここでは、モナドが持っている、計算をシリアライズするという性質だけが利用（濫用）されているに過ぎない。

Repa配列の型

というわけで、たとえばDoubleを即値で格納した二次元配列ならArray U (Z:.Int:.Int) Doubleという型になる。もちろんShape型のところは型シノニムを使ってArray U DIM2 Doubleと書いても差し支えない。

Repa配列を作成する

普通に作る　

これは別に難しくない。IArrayとかと似たようなものだが、インデックスがShape型クラスのインスタンスなのと、格納されている値がUとかBとかDとか色々あるので、それに応じた函数が使われる：

a :: Array U DIM2 Int 
a = fromListUnboxed (Z:.10:.5) [0..49]

みたいな感じ。配列の各次元の長さをShape型クラスのインスタンスの値で与える。

インデックスでアクセスするには

MapやIArrayでもお馴染みの(!)演算子を使うだけ。ただ、しつこいがインデックスの型がShape型クラスのインスタンスなことに注意。

a ! (Z:.9:.4)  ==>  49

始端が自由に決められるIArrayなどとは異なり、Repaでは 配列の始端はゼロ始まり なので、上記の配列にアクセスするときのインデックスはDIM2型の値(Z:.0:.0)から(Z:.9:.4)を使うことに注意。この範囲を超えてa ! (Z:.10:.5)とかすると、境界外でエラーになる(これはRepa配列の実体であるData.Vectorでの境界エラーになる)。というわけで、実体がVector値なので、リストからではなくVector値からももちろんRepa配列を生成できる（使う函数の名前はたとえば上記のものから"List"が取れたfromUnboxedとかになる）。

外部データから読み込む

配列データをファイルから読み込んだり、外部ポインタ(Foreign.Ptr)の先から読み取ることもできる。その辺のユーテリティはrepa-ioとかにある。多くの画像形式に対応したものはrepa-devilsにある。ここでは省略。

配列を操作する

mapとかfoldとかそういうの

mapは直感的に使える。注意すべきは返ってくる値がD型の配列になること。D型の配列は（中身がインデックスから値への函数なので）Show型クラスのインスタンスではないが、インデックスでアクセスすれば、サンクが計算されて値が取り出される。既出だが、計算を実行して配列全体を即値型にしたければComputeSかComputePを使う。一番内側の次元について畳み込んで次元をひとつ下げる函数としてfoldSとfoldPが使える。後者が結果をモナドに包んで返してくるのも一緒。

形状も変えるような一般的配列変換

これには｀traverse`を使う。これはやたらと複雑なので、チュートリアルで利用例を直接に参照すべき。

traverse （変換元の配列）（インデックス形状の変換函数）（要素値の変換函数）

インデックスについては次元数や要素数など形状を変えられる。値を変換する部分は、IArrayのixmapと似ていて、変換先の配列のインデックス値を取って、変換元の配列要素などを参照しつつ、変換先の配列の要素値になるべき値を返してくる（つまり全体的にcontravariantなmapを一般化したものとして作用する）。

算術演算

配列の要素ごとの算術演算のために、(+^),(*^)とかが用意されている、(*^)は`要素の乗算に過ぎないので、行列の積では「ない」ことに注意すべき。

計算の最適化とか

色々あるみたいだが、文系情弱には縁がないので省略。