Numeric Haskell: A Repa Tutorial

From HaskellWiki

Repa is a Haskell library for
high performance, regular, multi-dimensional parallel arrays. All
numeric data is stored unboxed and functions written with the Repa
combinators are automatically parallel (provided you supply "+RTS -N" on
the command line when running the program).

This document provides a tutorial on array programming in Haskell using the repa package.

1.2 Index types and shapes

Before we can get started manipulating arrays, we need a grasp of repa's notion of array shape. Much like the classic 'array' library in Haskell, repa-based arrays are parameterized via a type which determines the dimension of the array, and the type of its index. However, while classic arrays take tuples to represent multiple dimensions, Repa arrays use a richer type language for describing multi-dimensional array indices and shapes (technically, a heterogeneous snoc list).

Shape types are built somewhat like lists. The constructor Z corresponds
to a rank zero shape, and is used to mark the end of the list. The :. constructor adds additional dimensions to the shape. So, for example, the shape:

is the type of a two-dimensional array of doubles, indexed via `Int` keys, while

Array Z Double

is a zero-dimension object (i.e. a point) holding a Double.

Many operations over arrays are polymorphic in the shape / dimension
component. Others require operating on the shape itself, rather than
the array. A typeclass, Shape, lets us operate uniformally
over arrays with different shape.

1.2.1 Building shapes

To build values of `shape` type, we can use the Z and :. constructors. Open the ghci and import Repa:

Prelude>:m +Data.Array.Repa
Repa> Z -- the zero-dimension
Z

For arrays of non-zero dimension, we must give a size. Note: a common error is to leave off the type of the size.

Repa>:t Z :. 10
Z :. 10::Numhead=> Z :. head

leading to annoying type errors about unresolved instances, such as:

No instance for (Shape (Z :. head0))

To select the correct instance, you will need to annotate the size literals with their concrete type:

Repa>:t Z :. (10::Int)
Z :. (10::Int):: Z :. Int

is the shape of 1D arrays of length 10, indexed via Ints.

Given an array, you can always find its shape by calling extent.

Additional convenience types for selecting particular parts of a shape are also provided (All, Any, Slice etc.) which are covered later in the tutorial.

Repa arrays are instances of the Num. This means that
operations on numerical elements are lifted automagically onto arrays of
such elements. For example, (+) on two double values corresponds to
element-wise addition, (+), of the two arrays of doubles:

The swap function reorders the index space of the array.
To do this, we extract the current shape of the array, and write a function
that maps the index space from the old array to the new array. That index space function
is then passed to backpermute which actually constructs the new
array from the old one.

backpermute generates a new array from an old, when given the new shape, and a
function that translates between the index space of each array (i.e. a shape
transformer).

1.7.2 Example: matrix-matrix multiplication

A more advanced example from the Repa paper: matrix-matrix multiplication: the result of
matrix multiplication is a matrix whose elements are found by
multiplying the elements of each row from the first matrix by the
associated elements of the same column from the second matrix and
summing the result.

if and

then

So we take two, 2D arrays and generate a new array, using our transpose
function from earlier:

The idea is to expand both 2D argument arrays into 3D arrays by
replicating them across a new axis. The front face of the cuboid that
results represents the array a, which we replicate as often
as b has columns (colsB), producing
aRepl.

The top face represents t (the transposed b), which we
replicate as often as a has rows (rowsA), producing
bRepl,. The two replicated arrays have the same extent,
which corresponds to the index space of matrix multiplication

Optimized implementations of this function are available in the
repa-algorithms package.

1.7.3 Example: parallel image desaturation

To convert an image from color to greyscale, we can use the luminosity method to averge RGB pixels into a common grey value, where the average is weighted for human perception of green