repaの基本的な配列表現

　repaでは，並列配列としてData.Array.Repa.ReprモジュールでArray型が提供されています。repaのArray型は，第19回で説明した配列と名前は同じですが，全く異なる形で定義された別の型です。

-- | Arrays with a representation tag, shape, and element type.
data family Array r sh e

　repaが提供するArray型は，配列のデータを保持する内部表現を表すタグ，配列の要素を参照するために外部表現として使われる「行列の形状（shape）」，配列に格納される要素，の三つの型から構成されています。

　また，repaのArray型は，通常のデータ型としてではなく，第36回のコラムで説明した型族を使って定義されています。タグとして用意された型を型族の型変数rに渡すことで，Array型は対応した内部表現を持つ並列配列になります。

　例えば，Array型にU型のタグを渡すことで，非ボックス化Vectorを内部表現として使う配列になります。

-- | Unboxed arrays are represented as unboxed vectors.
--
--   The implementation uses @Data.Vector.Unboxed@ which is based on type
--   families and picks an efficient, specialised representation for every
--   element type. In particular, unboxed vectors of pairs are represented
--   as pairs of unboxed vectors.
--   This is the most efficient representation for numerical data.
--
data U
data instance U.Unbox e => Array U sh e
        = AUnboxed sh !(U.Vector e)

　Array型にD型のタグを渡すと，実データを持つ配列ではなく，「配列に対して計算処理を行う関数」を内部表現として使う，遅延評価の配列になります。

-- | Delayed arrays are represented as functions from the index to element value.
--
--   Every time you index into a delayed array the element at that position 
--   is recomputed.
data D
data instance Array D sh e
        = ADelayed  
                sh 
                (sh -> e)

　次に，Array型の2番目の引数として渡される，行列の形状について説明しましょう。行列の形状は，Shapeクラスのインスタンスとして定義されたZ型と「:. Int」型を使った「Z :. Int」，「Z :. Int :. ... Int」という型で表されます。

Prelude Data.Array.Repa> :i Shape 
class Eq sh => Shape sh where
  rank :: sh -> Int
  zeroDim :: sh
  unitDim :: sh
  intersectDim :: sh -> sh -> sh
  addDim :: sh -> sh -> sh
  size :: sh -> Int
  sizeIsValid :: sh -> Bool
  toIndex :: sh -> sh -> Int
  fromIndex :: sh -> Int -> sh
  inShapeRange :: sh -> sh -> sh -> Bool
  listOfShape :: sh -> [Int]
  shapeOfList :: [Int] -> sh
  deepSeq :: sh -> a -> a
        -- Defined in `Data.Array.Repa.Shape'
instance Shape Z -- Defined in `Data.Array.Repa.Index'
instance Shape sh => Shape (sh :. Int)
  -- Defined in `Data.Array.Repa.Index'

-- | An index of dimension zero
data Z = Z
       deriving (Show, Eq, Ord)

-- | Our index type, used for both shapes and indices.
infixl 3 :.
data tail :. head
        = !tail :. !head
        deriving (Show, Eq, Ord)

　例えば，4個の要素を持つ1次元の行列の形状は「Z :. (4 ::Int)」，2×3の2次元の行列の形状は「Z :. (2::Int) :. (3::Int)」，2×3×4の3次元の行列の形状は「Z :. (2::Int) :. (3::Int) :. (4::Int)」になります。

　Shapeクラスのrankメソッドは渡した形状の行列の次元を返す関数，sizeメソッドは渡した形状の行列に格納可能な要素数を返す関数です。これらの関数を適用してみましょう。

Prelude Data.Array.Repa> rank (Z :. (4::Int))
1
Prelude Data.Array.Repa> size (Z :. (4::Int))
4
Prelude Data.Array.Repa> rank (Z :. (2::Int) :. (3::Int))
2
Prelude Data.Array.Repa> size (Z :. (2::Int) :. (3::Int))
6
Prelude Data.Array.Repa> rank (Z :. (2::Int) :. (3::Int) :. (4::Int))
3
Prelude Data.Array.Repa> size (Z :. (2::Int) :. (3::Int) :. (4::Int))
24

　期待通りの結果ですね。

　実際には，いちいち「Z :. Int :. ... Int」という型や「Z :. (2::Int) :. ... (4::Int)」という値で行列の形状を表すのは煩わしいでしょう。そこで，行列の形状を表すのに便利な，いくつかの別名や補助関数が用意されています。

　プログラマーが行列を扱う場合に示したい情報は，「Z :. Int :. ... :. Int」のような構造そのものではありません。重要なのは「Z :. Int :. ... :. Int」という構造が，何次元の配列を表現しているかという情報です。

　そこでrepaでは，0次元から5次元までの行列の形状を表す型に対し，0から5に対応する数字nを持つDIMnという別名を用意しています。

-- Common dimensions
type DIM0       = Z
type DIM1       = DIM0 :. Int
type DIM2       = DIM1 :. Int
type DIM3       = DIM2 :. Int
type DIM4       = DIM3 :. Int
type DIM5       = DIM4 :. Int

　また，1次元から5次元までの行列の形状を表す値を簡単に作成できるよう，作成したい形状の行列の次元数Nを持つixN関数が補助関数として提供されています。

-- | Helper for index construction.
--
--   Use this instead of explicit constructors like @(Z :. (x :: Int))@.
--   The this is sometimes needed to ensure that 'x' is constrained to 
--   be in @Int@.
ix1 :: Int -> DIM1
ix1 x = Z :. x
{-# INLINE ix1 #-}

ix2 :: Int -> Int -> DIM2
ix2 y x = Z :. y :. x
{-# INLINE ix2 #-}

ix3 :: Int -> Int -> Int -> DIM3
ix3 z y x = Z :. z :. y :. x
{-# INLINE ix3 #-}

ix4 :: Int -> Int -> Int -> Int -> DIM4
ix4 a z y x = Z :. a :. z :. y :. x
{-# INLINE ix4 #-}

ix5 :: Int -> Int -> Int -> Int -> Int -> DIM5
ix5 b a z y x = Z :. b :. a :. z :. y :. x
{-# INLINE ix5 #-}

　では，実際にArray型にU型のタグを渡して配列を作ってみましょう。fromListUnboxed関数に行列の形状を渡すことで，Array型にU型のタグを渡した配列をリストから作ることができます。同様に，fromUnboxed関数に行列の形状を渡すことで，非ボックス化Vectorから作ることもできます。なお，「Array型にU型のタグを渡した配列」と記述するのは煩雑なので，以降は「U型の配列」と記述します。他のタグの場合も同様です。

Prelude Data.Array.Repa> :t fromListUnboxed  
fromListUnboxed
  :: (Shape sh, Data.Vector.Unboxed.Base.Unbox a) =>
     sh -> [a] -> Array U sh a
Prelude Data.Array.Repa> :t fromUnboxed  
fromUnboxed
  :: (Shape sh, Data.Vector.Unboxed.Base.Unbox e) =>
     sh -> Data.Vector.Unboxed.Base.Vector e -> Array U sh e

Prelude Data.Array.Repa> fromListUnboxed (Z :.4) [1..4] :: Array U DIM1 Double
AUnboxed (Z :. 4) (fromList [1.0,2.0,3.0,4.0])
Prelude Data.Array.Repa> fromListUnboxed (Z :.4) [1..4] :: Array U DIM1 Int
AUnboxed (Z :. 4) (fromList [1,2,3,4])
Prelude Data.Array.Repa> fromListUnboxed (Z :.(4::Int)) [1..4]
AUnboxed (Z :. 4) (fromList [1.0,2.0,3.0,4.0])
Prelude Data.Array.Repa> fromListUnboxed (Z :.(6::Int)) ['a'..'f']
AUnboxed (Z :. 6) (fromList "abcdef")
Prelude Data.Array.Repa> fromListUnboxed (Z :. 2 :. 3) ['a'..'f'] :: Array U DIM2 Char
AUnboxed ((Z :. 2) :. 3) (fromList "abcdef")
Prelude Data.Array.Repa> fromListUnboxed (ix1 4) [1..4]
AUnboxed (Z :. 4) (fromList [1.0,2.0,3.0,4.0])
Prelude Data.Array.Repa> fromListUnboxed (ix2 2 2) [1..4]
AUnboxed ((Z :. 2) :. 2) (fromList [1.0,2.0,3.0,4.0])
Prelude Data.Array.Repa> fromListUnboxed (ix3 2 3 4) ['a'..'x']
AUnboxed (((Z :. 2) :. 3) :. 4) (fromList "abcdefghijklmnopqrstuvwx")

Prelude Data.Array.Repa Data.Vector.Unboxed> let x = fromList [1..4]
Prelude Data.Array.Repa Data.Vector.Unboxed> fromUnboxed (Z :. 4) x :: Array U DIM1 Double
AUnboxed (Z :. 4) (fromList [1.0,2.0,3.0,4.0])
Prelude Data.Array.Repa Data.Vector.Unboxed> fromUnboxed (Z :. (4::Int)) x
AUnboxed (Z :. 4) (fromList [1.0,2.0,3.0,4.0])
Prelude Data.Array.Repa Data.Vector.Unboxed> let y = fromList [1..4] :: Vector Int
Prelude Data.Array.Repa Data.Vector.Unboxed> fromUnboxed (Z :. 2 :. 2) x :: Array U DIM2 Int
AUnboxed ((Z :. 2) :. 2) (fromList [1,2,3,4])
Prelude Data.Array.Repa Data.Vector.Unboxed> fromUnboxed (Z :. (2::Int) :. (2::Int)) x
AUnboxed ((Z :. 2) :. 2) (fromList [1.0,2.0,3.0,4.0])
Prelude Data.Array.Repa Data.Vector.Unboxed> fromUnboxed (ix1 4) $ fromList [1..4]
AUnboxed (Z :. 4) (fromList [1.0,2.0,3.0,4.0])
Prelude Data.Array.Repa Data.Vector.Unboxed> fromUnboxed (ix2 2 2) $ fromList [1..4]
AUnboxed ((Z :. 2) :. 2) (fromList [1.0,2.0,3.0,4.0])
Prelude Data.Array.Repa Data.Vector.Unboxed> fromUnboxed (ix3 2 2 2) $ fromList [1..8]
AUnboxed (((Z :. 2) :. 2) :. 2) (fromList [1.0,2.0,3.0,4.0,5.0,6.0,7.0,8.0])
Prelude Data.Array.Repa Data.Vector.Unboxed> fromUnboxed (ix3 2 3 4) $ fromList ['a'..'x']
AUnboxed (((Z :. 2) :. 3) :. 4) (fromList "abcdefghijklmnopqrstuvwx")

　逆に，toList関数を使うことで，Array型を使って表現される任意の配列からリストを作成できます。

Prelude Data.Array.Repa> :t toList 
toList :: (Shape sh, Repr r e) => Array r sh e -> [e]
Prelude Data.Array.Repa> let xs = fromListUnboxed (Z :.4) [1..4] :: Array U DIM1 Double
Prelude Data.Array.Repa> toList xs
[1.0,2.0,3.0,4.0]
Prelude Data.Array.Repa> let ys = fromListUnboxed (Z :. 2 :. 2) [1..4] :: Array U DIM2 Int
Prelude Data.Array.Repa> toList ys
[1,2,3,4]
Prelude Data.Array.Repa> let zs = fromListUnboxed (ix3 2 3 4 ) ['a'..'x']
Prelude Data.Array.Repa> toList zs
"abcdefghijklmnopqrstuvwx"

　また，toUnboxed関数を使うことで，U型の配列から非ボックス化Vectorを作成できます。

-- | O(1). Unpack an unboxed vector from an array.
toUnboxed
        :: U.Unbox e
        => Array U sh e -> U.Vector e
toUnboxed (AUnboxed _ vec)
        = vec
{-# INLINE toUnboxed #-}

Prelude Data.Array.Repa Data.Vector.Unboxed> :t toUnboxed xs
toUnboxed xs :: Vector Double
Prelude Data.Array.Repa Data.Vector.Unboxed> toUnboxed xs
fromList [1.0,2.0,3.0,4.0]
Prelude Data.Array.Repa Data.Vector.Unboxed> toUnboxed ys
fromList [1,2,3,4]
Prelude Data.Array.Repa Data.Vector.Unboxed> toUnboxed zs
fromList "abcdefghijklmnopqrstuvwx"

　U型の配列を作成する際には，作成する配列に格納できる要素数と，実際に格納される要素数を等しくする必要がある点に注意してください。fromListUnboxed関数では，作成する配列に格納できる要素数の定義と実際に格納される要素数が異なる場合にはエラーが生じます。

Prelude Data.Array.Repa> fromListUnboxed (ix3 3 6 5) ['a'..'z']
*** Exception: Data.Array.Repa.Eval.Fill.fromList: provide array shape does not match list length
Prelude Data.Array.Repa> fromListUnboxed (ix2 3 4) [0..99]
*** Exception: Data.Array.Repa.Eval.Fill.fromList: provide array shape does not match list length

　一方，fromUnboxed関数では，配列に格納できる要素数の定義と，内部表現として使われる非ボックス化Vectorに実際に格納される要素数が異なっていてもエラーになりません。

Prelude Data.Array.Repa Data.Vector.Unboxed> fromUnboxed (ix2 2 2) $ fromList [1..27]
AUnboxed ((Z :. 2) :. 2) (fromList [1.0,2.0,3.0,4.0,5.0,6.0,7.0,8.0,9.0,10.0,11.0,12.0,13.
0,14.0,15.0,16.0,17.0,18.0,19.0,20.0,21.0,22.0,23.0,24.0,25.0,26.0,27.0])
Prelude Data.Array.Repa Data.Vector.Unboxed> fromUnboxed (ix2 3 4) $ enumFromN 1 2
AUnboxed ((Z :. 2) :. 2) (fromList [1.0,2.0])
Prelude Data.Array.Repa Data.Vector.Unboxed> fromUnboxed (ix3 4 3 2) $ enumFromN 2 10
AUnboxed (((Z :. 4) :. 3) :. 2) (fromList [2.0,3.0,4.0,5.0,6.0,7.0,8.0,9.0,10.0,11.0])

　もっとも，配列に格納できる要素数の定義よりも実際に格納される要素数が小さい場合には，配列の要素を参照する際にエラーが生じます。

Prelude Data.Array.Repa Data.Vector.Unboxed> let x = fromUnboxed (ix2 3 4) $ enumFromN 1 2
Prelude Data.Array.Repa Data.Vector.Unboxed> :m -Data.Vector.Unboxed
Prelude Data.Array.Repa> index x (ix2 2 2)
*** Exception: ./Data/Vector/Generic.hs:244 ((!)): index out of bounds (10,2)
Prelude Data.Array.Repa> x ! (ix2 2 3)
*** Exception: ./Data/Vector/Generic.hs:244 ((!)): index out of bounds (11,2)

　エラーの発生の遅延はバグの原因になるので，こうしたコードを書くのは避けるべきです。

　配列の要素を参照するまでエラーの発生が遅延されるのは，fromUnboxed関数では，渡した非ボックス化Vectorは単純にU型の配列としてラップされ，検査されないためです。

-- | O(1). Wrap an unboxed vector as an array.
fromUnboxed
        :: (Shape sh, U.Unbox e)
        => sh -> U.Vector e -> Array U sh e
fromUnboxed sh vec
        = AUnboxed sh vec
{-# INLINE fromUnboxed #-}

　非ボックス化Vectorを内部表現として使う配列を非ボックス化Vectorから作るのであれば，非ボックス化Vectorの中にあるすべての要素を取り出して新しい配列を作る必要はありません。配列に格納できる要素数が一致するかどうかの検査を省略すれば，U型の配列作成のコストを低くできます。この結果，配列の要素を参照するときまでエラーの発生が遅延されるのです。

　配列の要素の参照にはindex関数や!演算子を使います。配列を作成する際と同様に，index関数や!演算子を使った配列の参照でも，DIM*型（「Z :. Int :. ... :. Int」型）の値を使います。

Prelude Data.Array.Repa> :t index
index :: (Shape sh, Repr r e) => Array r sh e -> sh -> e
Prelude Data.Array.Repa> :t (!)
(!) :: (Shape sh, Repr r e) => Array r sh e -> sh -> e

　index関数や!演算子の型には，行列の形状を表すShapeとは別に，Reprという型クラスの文脈が付いています。Reprクラスは，U型の配列やD型の配列など，様々な内部表現を持つ配列に対する参照操作を抽象化するためのものです。

Prelude Data.Array.Repa.Repr.Unboxed Data.Array.Repa> :i Repr
class Repr r e where
  extent :: Shape sh => Array r sh e -> sh
  index :: Shape sh => Array r sh e -> sh -> e
  unsafeIndex :: Shape sh => Array r sh e -> sh -> e
  linearIndex :: Shape sh => Array r sh e -> Int -> e
  unsafeLinearIndex :: Shape sh => Array r sh e -> Int -> e
  deepSeqArray :: Shape sh => Array r sh e -> b -> b
        -- Defined in `repa-3.1.4.2:Data.Array.Repa.Base'
instance Repr D a -- Defined in `Data.Array.Repa.Repr.Delayed'
instance Unbox a => Repr U a
  -- Defined in `Data.Array.Repa.Repr.Unboxed'

　index関数はReprクラスのメソッド，!演算子はindex関数の別名として定義されています。

-- | O(1). Alias for `index`
(!) :: (Repr r e, Shape sh) => Array r sh e -> sh -> e
(!) = index

　このような形で定義されたindex関数や!演算子を使うことで，U型やD型といった「Reprクラスのインスタンスとして定義された型タグ」を使って定義された，様々な内部表現を持つ配列に対して参照操作を行えます。

　なお，配列の参照操作を使う際には注意点があります。配列の作成と配列の参照のいずれの場合でも，Shapeクラスのインスタンスとして定義された型を使いますが，配列の作成のために渡す値と，配列の参照操作のために渡す値には違いがあるということです。配列の作成時には表現したい行列の大きさそのものを渡すのに対し，配列の参照では0起原の値を使います。例えば，2次元のn×mという行列の場合，配列の作成には「ix2 n m」を使いますが，配列の要素の参照には「ix2 0 0」から「ix2 (n-1) (m-1)」までしか利用できません。このため，「ix2 n m」という値で配列を参照しようとすると，範囲外を参照しようとしたというエラーになってしまいます。

Prelude Data.Array.Repa> index (fromListUnboxed (ix2 4 5) [1..20]) (ix2 4 5)
*** Exception: ./Data/Vector/Generic.hs:244 ((!)): index out of bounds (25,20)

　では，U型の配列の要素を参照してみましょう。

Prelude Data.Array.Repa> let x = fromListUnboxed (ix2 3 4) [1..12]
Prelude Data.Array.Repa> index x (ix2 2 3)
12.0
Prelude Data.Array.Repa> x!(ix2 1 2)
7.0
Prelude Data.Array.Repa> let y = fromListUnboxed (ix3 3 4 5) [1..60]
Prelude Data.Array.Repa> index y (ix3 1 3 2)
38.0
Prelude Data.Array.Repa> y!(ix3 2 3 2)
58.0

　この例では配列の範囲内を指しているので，配列に格納した要素がきちんと返ってきています。

　次にD型の配列の要素を参照してみましょう。delay関数を使うことで，U型の配列からD型の配列を作成できます。

-- | O(1). Delay an array.
--   This wraps the internal representation to be a function from
--   indices to elements, so consumers don't need to worry about
--   what the previous representation was.
--
delay   :: (Shape sh, Repr r e)
        => Array r sh e -> Array D sh e
delay arr = ADelayed (extent arr) (unsafeIndex arr)
{-# INLINE delay #-}

　delay関数を使ってD型の配列を作り，その値を参照してみましょう。

Prelude Data.Array.Repa> let x' = delay x
Prelude Data.Array.Repa> index x' (ix2 2 3)
12.0
Prelude Data.Array.Repa> x'!(ix2 1 2)
7.0
Prelude Data.Array.Repa> let y' = delay y
Prelude Data.Array.Repa> index y' (ix3 1 3 2)
38.0
Prelude Data.Array.Repa> y'!(ix3 2 3 2)
58.0

　U型の値を参照したときと同じ値が返ってきているのがわかりますね。delay関数ではU型の配列をD型の配列に変換しただけで，その内部の要素に変更は加えていません。したがって，これは期待通りの振る舞いです。

　delay関数でD型の配列に変換できるのは，U型の配列だけではありません。delay関数はReprクラスのメソッドを使って定義されています。このためdelay関数は，特定の内部表現を持つ配列からだけではなく，「Reprクラスのインスタンスとして定義された型タグ」を使って定義された様々な内部表現を持つ任意の配列から，D型の配列を作成できます。例えば，D型の配列からさらにD型の配列を作ることもできます。

Prelude Data.Array.Repa> let y'' = delay y'
Prelude Data.Array.Repa> index y'' (ix3 1 3 2)
38.0
Prelude Data.Array.Repa> y''!(ix3 2 3 2)
58.0

　内部の要素を変更しているわけではないので，D型の配列から作ったD型の配列もやはり同じ値を返します。