
\section{The main new idea: Type Indexing}
\label{sec:type-indexing}

We are now ready to explain our main technical innovation.
In Repa~version~3 we define arrays as a \emph{data family} \cite{Chakravarty:AssocTypes}:
%
\begin{code}
    data family Array rep sh e
\end{code}
%
An @Array@ represents a partial function from indices of type @sh@ to elements of type @e@.  The array is defined on a range of indices, from zero to a maximum called the \emph{extent} of the array. In this family, @rep@ is a \emph{type index} that specifies the representation of the array, while @sh@ is the shape, and @e@ is the element type as before. Figure~\ref{figure:Repa3Array} gives two particular instances of @Array@, where @D@ is for delayed and @U@ is for (unboxed) manifest arrays respectively.

\eject
% -----------------------------------------------------------------------------
\begin{figure}
\begin{small}
\begin{code}
data family Array rep sh e
data instance Array D sh e = ADelayed sh (sh -> e)
data instance Array U sh e = AUnboxed sh (Vector e)
...etc...

-- The type indices are declared as nullary types
data D   -- Delayed
data U   -- Manifest, unboxed
...etc...
\end{code}
\end{small}
\caption{Towards the Repa 3 Array Representation}
\label{figure:Repa3Array}
\end{figure}

We will give more detail shortly, but we can already see how type indexing addresses the problems of \S\ref{section:problem}:

\begin{itemize}
\item We can give a more informative type to @force@:
\begin{small}
\begin{code}
 force :: Shape sh => Array D sh e -> Array U sh e
\end{code}
\end{small}
Unlike \S\ref{section:force}, this type statically specifies that the input is delayed and the output is manifest. We cannot accidentally force a manifest array, or forget to force a delayed array in a situation where a manifest one is needed.

\item A function like @(f :: Array U sh e -> ...)@ receives an argument that can only be built with @AUnboxed@. There is no redundant tag-testing to check that the array has the representation that the programer already knows it has (\S\ref{section:rep-tests}).

\item When there are exactly (say) three regions in an array, we can use a \emph{type-level list} to drive code specialisation, avoiding the ad-hoc approach of \S\ref{section:Explosion}.  Details in \S\ref{section:Partitioned}.
\end{itemize}
Better still, type indexing scales up to allow a variety of different array representations, each with a different storage/performance cost model. Indeed, Repa~3 has no fewer than ten such representations:
\begin{itemize}
\item @D@ -- Delayed arrays (delayed)            \S\ref{section:Source}
\item @C@ -- Cursored arrays (delayed)           \S\ref{section:Cursored}
\item @U@ -- Adaptive unboxed vectors (manifest) \S\ref{section:Source}
\item @V@ -- Boxed vectors (manifest)            \S\ref{section:ReprPoly}
\item @B@ -- Strict byte arrays (manifest)       \S\ref{section:ReprPoly}
\item @F@ -- Foreign memory buffers (manifest)   \S\ref{section:ReprPoly}
\item @P@ -- Partitioned arrays (meta)           \S\ref{section:Partitioned}
\item @S@ -- Smallness hints (meta)              \S\ref{section:Smallness}
\item @I@ -- Interleave hints (meta)             \S\ref{section:Interleaved}
\item @X@ -- Undefined arrays (meta)             \S\ref{section:Undefined}
\end{itemize}
We can think of the type indices being generated by this kind declaration:
\begin{small}
\begin{code}
  kind RepIndex = D | C | U | V | B | X
                | P RepIndex RepIndex
                | S RepIndex | I RepIndex
\end{code}
\end{small}
With this declaration, @Array (P U X) sh e@ is a valid array type. GHC's recent @DataKinds@ extension supports exactly this form of declaration. However, using data kinds would make the index kind \emph{closed}, preventing users from adding new representation indices. Instead, we proceed as in Figure~\ref{figure:Repa3Array}, declaring type indices (such as @D@ and @U@) as fresh uninhabited types. An open and extensible set of array representations enables integration with other array libraries as the integration with the Accelerate GPGPU library shows.\footnote{\url{https://github.com/AccelerateHS/accelerate-io}}

Of these, the @D@ and @C@ indices specify \emph{delayed} array representations, meaning the array is expressed as a function from (value) indices to array elements. In this paper we refer to cursored arrays as being ``delayed'' as well due to the nature of the representation.

The @U@, @V@, @B@ and @F@ indices specify \emph{manifest} representations, meaning real data in memory. Supporting multiple manifest representations makes it easier to interface with third-party array libraries, such as @bytestring@. The Foreign (@F@) representation allows us to compute array elements and store them directly into foreign memory buffers, perhaps provided by the operating system. This eliminates copying between the GHC heap and foreign memory that would otherwise be necessary.

Finally, the @P@, @S@, @I@ and @X@ indices specify \emph{meta} representations. They combine other array types or encode information that does not directly define array elements. The \emph{partitioned} (@P@) and \emph{undefined} (@X@) representations together provide the partitioned arrays from \S\ref{section:Partitioned}. The \emph{smallness hint} (@S@) ensures that an array is evaluated sequentially, and the \emph{interleave hint} (@I@) manage unbalanced workloads by making each thread compute alternate array elements.


% -----------------------------------------------------------------------------
\subsection{Representation-polymorphic Operations}
\label{section:Source}
\begin{figure}
\begin{small}
\begin{code}
-- The Source class ----------------------------------
class Source r e where
 data Array r sh e
 extent       :: Shape sh => Array r sh e -> sh
 index        :: Shape sh => Array r sh e -> sh  -> e
 linearIndex  :: Shape sh => Array r sh e -> Int -> e

instance Source D e where
 data Array D sh e = ADelayed !sh (sh -> e)
 extent       (ADelayed sh _)    = sh
 index        (ADelayed _  f) ix = f ix
 linearIndex  (ADelayed sh f) ix = f (fromIndex sh ix)

instance Vector.Unbox e => Source U e where 
 data Array U sh e = AUnboxed !sh !(Vector e)
 ...

-- The Target class ----------------------------------
class Target rt e where
 data MVec rt e
 newMVec          :: Int -> IO (MVec rt e)
 unsafeWriteMVec  :: MVec rt e -> Int -> e -> IO ()
 unsafeFreezeMVec :: sh  -> MVec rt e 
                         -> IO (Array rt sh e)

instance Vector.Unbox e => Target U e where
 data MVec U e = UMVec (IOVector e)
 ...

-- The Load Class ------------------------------------
class (Shape sh, Source rs e) => Load rs sh e where
 loadP, loadS :: Target rt e
              => Array rs sh e -> MVec rt e -> IO ()
\end{code}
\end{small}
\caption{Repa 3 Array Representation}
\label{figure:Source}
\label{figure:Target}
\end{figure}
If we know the type index, we know the array representation, but what if we don't? Fundamental operations, such as array indexing, ought to be polymorphic in the representation. Since the representation of arrays varies with the type index, all polymorphic operations must involve a corresponding type. The canonical approach is to make @Array@ an \emph{associated type} of a class \cite{Chakravarty:AssocTypes}. The resulting declarations are in Figure~\ref{figure:Source}, which replaces Figure~\ref{figure:Repa3Array}.

We see that @Array@ is now an associated type of @Source@.  For each type index (@D@, @U@, @V@, and so on) we define an instance of @Source@, and each such instance gives a @data@ instance declaration for @Array@. The methods of @Source@ allow us to perform representation-polymorphic operations on @Array@s. The @extent@ function takes the shape of an array, and @index@ provides shape-polymorphic indexing. The @linearIndex@ function accesses the underlying flat representation.

The @Source@ class must be parameterised over the representation index @r@ as well as the array element type @e@, because certain representations restrict the element type. In particular, to store array elements unboxed (type index @U@) we need to know (a) the width of the unboxed elements, and (b) the data constructor to use when boxing them. This information is encapsulated in the @Unbox@ class, defined by the standard @Vector@ library, and used in the instance declaration for @Source U e@ in Figure~\ref{figure:Source}.

Note that the @Source@ class contains operations that \emph{read} or \emph{consume} an array only. It does not offer operations that \emph{construct} an array. We keep array-construction methods separate because, in general, they depend on both the source and result representations, and thus require \emph{two} type indices. We discuss this next.


% -----------------------------------------------------------------------------
\subsection{Parallel Computation and Array Filling}
\label{section:ParallelComputation}

Repa represents manifest arrays by real data in memory. It constructs a manifest array by first allocating a new array, performing parallel writes to initialise the elements, and then freezing the result into an immutable version. These three operations\footnote{The latter two operations have names starting with @unsafe@ because @unsafWriteMVec@ does no bounds checking, and @unsafeFreezeMVec@ does not check that further writes do not take place. This is an internal interface that is not normally used by client programmers.} are bundled into the @Target@ class (Figure~\ref{figure:Target}). The @MVec@ associated type specifies the underlying type of one-dimensional mutable vectors. Delayed and meta representations do not correspond to real data in memory, thus cannot be written to and are not instances of @Target@.

The @Load@ class, also shown in Figure~\ref{figure:Target} forms the bridge between @Source@ and @Target@. The @loadP@ function of @Load@ takes an immutable source array of type @Array rs sh e@, and a mutable destination vector of type @MVec rt e@. Instances of @loadP@ will fork several threads to concurrently read elements from the source array, and write them to the destination. The @loadS@ function performs the same operation sequentially, which we discuss further in \S\ref{section:Smallness}. 

With @Load@ and @Target@, we can write the generic parallel array computation function, @computeP@, taking the role of @force@ (\S\ref{section:force}):
%
\begin{small}
\begin{code}
   computeP :: (Load rs sh e, Target rt e)
            => Array rs sh e -> Array rt sh e
   computeP arr1
    = unsafePerformIO
    $ do mvec2 <- newMVec (size $ extent arr1) 
         loadP arr1 mvec2
         unsafeFreezeMVec (extent arr1) mvec2
\end{code}
\end{small}
%
In Repa~3 we use the name @computeP@ instead of @force@, because the provided @Load@ instances only allow delayed and meta representations to be used as the source. With these representations, @loadP@ runs the delayed computations in parallel. To copy data between \emph{manifest} representations we provide a separate function @copyP@ which delays the source array before applying @computeP@. In Repa 1 \& 2, applying @force@ to an already @Manifest@ array was a no-op. This behaviour turned out to be unnecessary, and needlessly confusing for client programmers.

Finally, we keep the @Load@ and @Source@ classes separate because, for some array representations (@rs@), we want to provide a different @loadP@ instance for each shape (@sh@). Specifically, the @loadP@ function for @DIM2@ cursored arrays uses a column-based traversal order, which we discuss further in \S\ref{section:Interleaved}. 


% -----------------------------------------------------------------------------
\subsection{Bulk Operations}

\begin{figure}
\begin{small}
\begin{code}
delay :: (Shape sh, Source r e)
      => Array r sh e -> Array D sh e
delay arr = ADelayed (extent arr) (index arr)

map   :: (Shape sh, Source r a)
      => (a -> b) -> Array r sh a -> Array D sh b
map f arr = case delay arr of
               ADelayed sh g -> ADelayed sh (f . g)

zipWith :: (Shape sh, Source r1 a, Source r2 b)
        => (a -> b -> c) 
        -> Array r1 sh a -> Array r2 sh b -> Array D sh c
zipWith f arr1 arr2
 = let  ADelayed sh1 f1 = delay arr1
        ADelayed _   f2 = delay arr2
        get ix = f (f1 ix) (f2 ix)
   in   ADelayed sh1 get
\end{code}
\end{small}
\caption{Bulk operations}
\label{figure:BulkOperations}
\label{figure:map}
\end{figure}

The definitions of @map@ and @zipWith@ using our new setup are shown in Figure~\ref{figure:BulkOperations}. While the bodies of these functions are almost identical to the Repa 1 versions from Figure~\ref{figure:Repa1}, their types now express that the result has the Delayed @(D)@ representation. 

Redoing the @doubleZip@ example from \S\ref{section:doubleZip} yields:
%
\begin{small}
\begin{code}
 doubleZip1 :: Array U DIM2 Int -> Array U DIM2 Int 
            -> Array D DIM2 Int
 doubleZip1 arr1 arr2
  = map (* 2) $ zipWith (+) arr1 arr2
\end{code}
\end{small}
%
Here we have given a type signature that explicitly constrains the input arrays to have Unboxed @(U)@ representation. The result array must have Delayed (@D@) representation, corresponding to the signature of @map@. 

When we apply @computeP@ the result array must have a manifest representation, as only manifest representations can be made instances of @Target@. For example, we could constrain the result to have Unboxed (@U@) representation as well:
\par
\begin{small}
\begin{code}
  doubleZip2 :: Array U DIM2 Int -> Array U DIM2 Int 
             -> Array U DIM2 Int
  doubleZip2 arr1 arr2
   = computeP $ map (* 2) $ zipWith (+) arr1 arr2
\end{code}
\end{small}
%
There is no need to provide explicit patterns such as @Manifest{}@ for the parameter variables as we did in \S\ref{section:Introduction}, because the type index controls the representation directly. Alternately, if we also leave the array shape polymorphic, then the most general type we could assign is the following:
\par
\begin{small}
\begin{code}
  doubleZip3 :: ( Source r1 e, Source r2 e, Target r3 e
                , Shape sh, Num e)
             => Array r1 sh e -> Array r2 sh e
             -> Array r3 sh e
  doubleZip3 arr1 arr2
   = computeP $ map (* 2) $ zipWith (+) arr1 arr2
\end{code}
\end{small}
%
We now return to our stated goal of helping client programmers write fast programs. The fact that an array's representation is determined by its type index provides an easy-to-follow rule specifying when to attach @INLINE@ pragmas to client-written functions: 
\begin{quote}
If the signature contains any @D@ tags, or type class dictionaries, then you must @INLINE@ it to get fast code, otherwise not.
\end{quote}
Inlining forces the user-defined element functions in delayed arrays to be fused with applications of @computeP@ present in callers. It also ensures that the instance functions present in type class dictionaries are inlined and fused appropriately. For @doubleZip2@, we do not need to inline it into its caller to get fast code, because the only code fusion that takes place happens within the function itself.

