04-frontend.tex

%% We use `subfiles' package
\documentclass[preamble.tex]{subfiles}
\begin{document}

\clearpage

\chapter{Delaying array computations}
\label{ch:frontend}

This chapter describes \LiveFusion\ -- a new array language that I designed. It offers a library of fusible array combinators and is able to exploit more fusion opportunities than equational fusion frameworks\ieqf (see Appendix \ref{sec:Equational-Fusion-Systems})

In its essence \LiveFusion takes the approach of \*deeply embedded domain specific languages* (EDSLs).\iedsl


\section{On embedded domain-specific languages}
\iedsl

The term \term{domain-specific language} (DSL) describes a computer language designed for solving problems in a particular application domain.

Examples of DSLs include \name{HTML} for webpages markup \cite{HTML}, \name{SQL} for relational database queries \cite{SQL} and \name{Verilog} for hardware description \cite{Verilog}.

In contrast, \Haskell \cite{Pey03} is a \*general purpose language* and is more widely applicable across domains.

On the other hand, an \term{embedded domain-specific language} (EDSL) is a domain-specific language which is designed to be used from within another language, called \*host* language. It is implemented as a library for the host language and can reuse its syntax, compiler and the runtime system.

The added benefit of the approach is that the EDSL code can interact with the rest of the code in the host language. This way the programmers do not have to leave the familiar environment and can write EDSL programs which are seamlessly integrated with their programs.

Two types of embedding exist:
\begin{description}
\item \term{Shallow embedding}

Shallowly embedded DSLs directly translate their constructs to expressions in the host language. They offer a fixed interpretation of their operations.

\name{Parsec}\footnote{parsec: Monadic parser combinators - http://hackage.haskell.org/package/parsec} library of parsing combinators is an example of a shallow EDSL for \Haskell.

\item \term{Deep embedding} 

The principle idea behind \term{deeply embedded DSLs} is that all operations of such an EDSL construct a tree of operations, or an \term{abstract syntax tree} (AST).\iast This new AST is a data structure in the host language. It is separate from the host language's AST and is rather presented as a tree-like structure in the host language's runtime.

The difference from the ASTs found in compilers is that it allows for analysis, optimisation and compilation (or interpretation) of the EDSL's AST at runtime. The complete process is taken care by the library written in the host language. This avoids the difficulties of writing a complete standalone optimising compiler.

\name{Accelerate} language for functional programming on GPUs \cite{CKL+11} is a vivid example of a deeply embedded DSL.
\end{description}

\LiveFusion EDSL falls under the category of deeply embedded domain-specific languages. In this chapter I describe the process of embedding array computations.

Figure~\ref{fig:Overview} presents the overview of the \LiveFusion language and the way it fits into the \name{Data Parallel Haskell} framework. While many of the concepts in the figure have not yet been discussed, it can serve as reference to the reader in future chapters.

\begin{figure}
  \includegraphics[width=\textwidth]{img/Overview.pdf}
  \caption{Overview of \DPH and \LiveFusion systems.}
  \label{fig:Overview}
\end{figure}



\section{\LiveFusion EDSL by example}

Figure \ref{fig:Flat-FilterMax} presents a program that can be written using \LiveFusion library as well as the corresponding data flow diagram.

\begin{figure}
\begin{subfigure}{0.5\textwidth}
\begin{hscode}
{-# LANGUAGE OverloadedLists #-}

filterMax xs 
  = let xs' = map (+1) xs
        pos = filter (>. 0) xs'
        mx  = fold max 0 xs'
    in  (pos :*: mx)

main = print
     $ filterMax [1,-2,5,0]

-- Output: (fromList [2,6,1], 6)
\end{hscode}
\end{subfigure}
%
\begin{subfigure}{.5\textwidth}%
\includegraphics[center,scale=\omniscale]{img/DFD-Flat-FilterMax}%
\end{subfigure}%

\caption{Example of a \LiveFusion program (left) and its data flow diagram (right).}
\label{fig:Flat-FilterMax}
\end{figure}

The function uses the familiar @map@, @filter@ and @fold@ combinators. Operators and functions of @Num@ type class (@+@, @-@, @abs@...) and number literals (0, 0.5, 1...) are supported directly. Versions of @Ord@ and @Eq@ operators are available postfixed with a dot (@>.@, @==.@ ...).

The example employs the explicit tuple constructor @(:*:)@, which brings together the branches computing @pos@ and @mx@ to be analysed for fusion. It addresses the problem of \*tupling* that prohibits fusion of independent terms in a pure context.

A vectorised version of this function will be presented in the Section~\ref{sec:FilterMax} when we look at the \QuickHull algorithm implementation in \DPH.

%\LiveFusion maps most of its combinators to AST nodes discussed next.%It is now time to show the AST node constructors these combinators correspond to. 



\section{\LiveFusion interface and AST}
\label{sec:AST}

\LiveFusion EDSL offers a library of flat\iflatcomb and segmented array combinators\isegcomb some of which are presented in Listing \ref{lst:LiveFusion-interface}.

\begin{hscode2}[%
    caption={\LiveFusion interface functions. Typeclass constraints on array elements have been omitted for brevity.},%
    label=lst:LiveFusion-interface,%
]
type Array  a = ArrayAST  a
type Scalar a = ScalarAST a

%\CommentPlusPlus{Flat array combinators}%
map      :: (Term a -> Term b) -> Array a -> Array b
filter   :: (Term a -> Term Bool) -> Array a -> Array a
zipWith  :: (Term a -> Term b -> Term c) -> Array a -> Array b -> Array c
zip      :: Array a -> Array b -> Array (a,b)
fold     :: (Term a -> Term a -> Term a) -> Term a -> Array a -> Scalar a
scan     :: (Term a -> Term a -> Term a) -> Term a -> Array a -> Array a
replicate:: Term Int -> Term a -> Array a
bpermute :: Array a      -- Source array
         -> Array Int    -- Indices to pick elements from
         -> Array a
packByTag:: Term Bool -> Array Bool -> Array a -> Array a

%\CommentPlusPlus{Segmented array combinators}%
scan_s   :: (Term a -> Term a -> Term a) -> Term a
         -> Array Int -> Array a -> Array a
fold_s   :: (Term a -> Term a -> Term a) -> Term a
         -> Array Int -> Array a -> Array a
replicate_s :: Term Int -> Array Int -> Array a -> Array a
         % %
\end{hscode2}


\begin{hscode2}[%
    caption={\LiveFusion AST (partial).},%
    label=lst:LiveFusion-AST,%
]
type ArrayAST  a = AST (Vector a)
type ScalarAST a = AST a

data AST t where
  %\CommentPlusPlus[21em]{Flat combinators}% 
  Map         :: (Elt a, Elt b)
              => (Term a -> Term b)
              -> ArrayAST a
              -> ArrayAST b

  Filter      :: Elt a
              => (Term a -> Term Bool)
              -> ArrayAST a
              -> ArrayAST a

  Fold        :: Elt a
              => (Term a -> Term a -> Term a)
              -> Term a
              -> ArrayAST a
              -> ScalarAST a

  Replicate   :: Elt a
              => Term Int
              -> Term a
              -> ArrayAST a
  %\CommentPlusPlus[21em]{Segmented combinators}%
  Fold_s      :: Elt a
              => (Term a -> Term a -> Term a)
              -> Term a
              -> ArrayAST Int
              -> ArrayAST a
              -> ArrayAST a

  Replicate_s :: Elt a
              => Term Int
              -> ArrayAST Int
              -> ArrayAST a
              -> ArrayAST a
  %\CommentPlusPlus[21em]{Manifest arrays}%
  Manifest    :: Elt a
              => Vector a
              -> ArrayAST a
  %\CommentPlusPlus[21em]{Explicit tupling}%
  (:*:)       :: (Typeable t1, Typeable t2)
              => AST t1
              -> AST t2
              -> AST (t1,t2)
           % %
\end{hscode2}

The fundamental concept by which \LiveFusion makes fusion possible is constructing an abstract syntax tree (AST)\iast of pending array operations at runtime and compiling that AST to efficient code when the result is required in the host program. Thus the compile time of \LiveFusion language maps to the runtime of its host language -- \Haskell.%By constructing an AST in a running program, the actual computations are \*delayed* to be performed at a later time.

The AST is the topmost layer of \LiveFusion system. Most of the library's user-facing functions are just constructing the nodes of the AST. Parts of \LiveFusion AST are presented in Listing~\ref{lst:LiveFusion-AST}. Is it easy to spot the correspondence between the user-facing functions from Listing~\ref{lst:LiveFusion-interface} and the constructors of the AST.

\LiveFusion AST is implemented as a \*Generalised Algebraic Data Type* (GADT) \cite{Jones:2006eh} in order to have more control over the types of individual constructors.

In principle, given an interpreting function such as @eval :: AST t -> t@, a value of type @t@ encoded in the @AST@ language can be computed. Notably, @t@ in @AST t@ is not necessarily an array. In fact not all of the array combinators return an array: in particular @fold@ and several others return a scalar value.


\subsection{Sharing recovery}
\label{sec:sharing-recovery}

In pure functional languages like \Haskell sharing of terms is \*not observable*. That is, given two terms it is not possible to tell whether or not they reside at the same memory address (as it is possible with pointers in a procedural language). Thus, if \LiveFusion AST constructed for the earlier @filterMax@ example (Figure~\ref{fig:Flat-FilterMax}) was traversed, @Map@ and @Manifest@ nodes would be encountered twice as in Figure~\ref{fig:Flat-FilterMax-Sharing} (left). Naively compiling such code would result in two computations of the result of @Map@.\footnote{Notably, this is an instance of fusion into multiple consumers identified as a major challenge for a fusion system to overcome in Section~\ref{sec:multiple-consumers}.}

Nonetheless, implementing support for sharing in an EDSL is possible. An overview of several approaches is given in \cite{ImplicitExplicitSharing}. There is a distinction between \*explicit* and \*implicit* sharing. When using explicit sharing, the sharing of terms is specified by the programmer using some explicit construct.

Moving on to \*implicit* sharing, there are techniques based of \*structural equality* of terms where two subtree of terms are considered shared if they are equal in structure and values. This approach is possible to implement without leaving the pure context (provided terms can be compared for equality). However, the experience of \Accelerate language authors \cite{CKL+11} shows that this approach suffers from severe complexity blowups for sufficiently large problems.

The second group of \*implicit* sharing approaches involves using techniques to escape the pure context to acquire pointers or otherwise stable references to a term. Sharing recovery using this approach grows linearly in the size of the tree. \LiveFusion incorporates sharing recovery of this type. The approach is called \name{Observable Sharing} and is described in \cite{ObservableSharing}.

\begin{figure}
\includegraphics[scale=0.7]{img/AST-Flat-FilterMax-Unshared}%
\includegraphics[scale=0.7]{img/AST-Flat-FilterMax-Shared}
\caption{\*Left:* Conceptual representation of \code{filterMax} AST in a pure context.\\\*Right:* Likely runtime representation of \code{filterMax} node graph. Desired result of sharing recovery.}
\label{fig:Flat-FilterMax-Sharing}
\end{figure}

More specifically, the approach is adapted to support a GADT AST \cite{GADTSharing}. An \*abstract semantic graph (ASG)* data type is specified (Listing~\ref{lst:ASG}) which has the same constructors as the original AST in addition to @VarG@ -- a representation of a shared variable in the graph. Then a @mapDeRef@ function is implemented that converts the AST to ASG, \*replacing every point of recursion with a* @VarG@ node. Sharing recovery of @filterMax@ results in the following graph, where recursive variable references have been replaced with unique integers:


\begin{hscode}
-- Nodes:
[
  (1,Wrap (BothG (VarG 2) (VarG 5))),
  (2,Wrap (FilterG <function> (VarG 3))),
  (3,Wrap (MapG <function> (VarG 4))),
  (4,Wrap (ManifestG (fromList [1,-2,5,0]))),
  (5,Wrap (FoldG <function> 0 (VarG 3))),
]
-- Entry node:
1
\end{hscode}


While the approaches to sharing recovery based on obtaining unique references to the term are generally non-portable and unsafe, it is guaranteed to be free of false positives. That is no erroneous sharing of terms that are not shared will occur. While false negatives are a possibility (i.e. unrecovered sharing), they are rare in practice.

\begin{hscode2}[%
  caption={Abstract semantic graph type for \LiveFusion AST and AST conversion function.},%
  label={lst:ASG}%
]
data ASG t s where
  MapG      :: (Elt a, Elt b)
            => (Term a -> Term b)
            -> ArrayASG a s
            -> ArrayASG b s
  ...

  VarG      :: Typeable t
            => s
            -> ASG t s


instance Typeable t => MuRef (AST t) where
  type DeRef (AST t) = WrappedASG
  mapDeRef ap t = Wrap <$> mapDeRef' ap t
    where
      mapDeRef' :: Applicative ap
                => (forall b. (MuRef b, WrappedASG ~ DeRef b) => b -> ap u)
                -> AST t
                -> ap (ASG t u)

      mapDeRef' ap (Map f arr)
        = MapG f
          <$> (VarG <$> ap arr)

      mapDeRef' ...
\end{hscode2}


\section{Parametrising higher-order combinators}
\label{sec:Scalar-language}
\iscalarlang

Another notable feature of the interface to note here is the type of functions that parametrise higher-order combinators:

\begin{hscode}
map    :: Elt a => (Term a -> Term b) -> Array a -> Array b
filter :: Elt a => (Term a -> Term Bool) -> Array a -> Array Bool
fold   :: Elt a => (Term a -> Term a -> Term a) -> Term a -> Array a
                                                          -> Scalar a
\end{hscode}

The corresponding \Haskell list functions or array libraries like @Vector@ or @Repa@ \cite{KCL+10} accept as arguments vanilla \Haskell functions like @a -> b@ or @a -> a -> a@ as long as the types match.

The question that should be asked, however, is why cannot \LiveFusion do the same?

%\LiveFusion AST is open to interpretation by several backends as discussed in Chapter \ref{ch:Code-Generation}.
The answer lies in the approach to evaluating the AST. In order to be efficient the AST needs to be compiled to highly optimised code. During compilation (discussed in Chapter \ref{ch:Code-Generation}) the AST is compiled \*at runtime* to one or more loops for which \Haskell code is generated.

If the parametrising functions were simple \Haskell functions, they would have been compiled to machine code \*before* the program is run. When the AST is later constructed at runtime they would only be available in the form of closures which can be called but not inspected\footnote{There is a non-portable way to explore closures in \Haskell, but it will not allow one to easily make use of the code compiled to binary.}.

Suppose the user writes @map (+1) xs@. The function @(+1)@ needs to be inlined into the compiled loop code to produce efficient code. If this does not happen, calling @(+1)@ for every element in the @xs@ array will result in the following:

\begin{enumerate}
\halfspacing
\item Boxing the element value by placing it in the heap
\item Calling the closure with the pointer to the boxed value
\item Jumping to the address of the @(+1)@ function referenced in the closure
\item The code will unbox the value, increment it and create a new value on the heap
\item The loop code will then need to unbox the result and write it into the result array
\end{enumerate}

This is a very round-about way of incrementing a value in a tight loop and will greatly affect the performance of the loop.

Unfortunately, at the time of generating code for AST, all these Haskell functions will have been compiled to machine code and are only available as closures. There is currently no way to inline them or reify them into the original \Haskell source code from which they were created.

I needed to find a way to generate efficient code for user specified functions at runtime. I needed a way to record the user-provided functions in a more high level way than machine code.

In summary the chosen representation for functions:
\begin{enumerate}
\halfspacing
\item Provides backend independent interface offering many common functions
\item Allows backend-specific implementation for each function
\item Hides implementation from the user by default
\item Allows user to compose functions in a way that looks native
\item Allows user to provide direct backend-specific implementation for new functions that cannon be composed from functions already provided
\end{enumerate}

The following section explains how this is achieved.


\subsection{Scalar Term language}

An expression of the type @Term a -> Term b@ is a term in the scalar expression language called \name{Term}, representing a function from type @a@ to type @b@. The AST for the \name{Term} language is given in Listing~\ref{lst:HOAS}.

\begin{hscode2}[%
    caption={Term language for HOAS representation.},%
    label=lst:HOAS,%
]
data Term t where
  -- Function or constant in backend-specific form
  Con :: Impl t -> Term t

  -- Lambda abstraction
  Lam :: (Term s -> Term t) -> Term (s -> t)

  -- Function application
  App :: Term (s -> t) -> Term s -> Term t

  -- for conversion to de Bruijn representation
  Tag :: Level -> Term t

type Level = Int
\end{hscode2}

The language allows the user to write functions in the higher-order abstract syntax (HOAS) form which looks native and in many cases reuses the familiar operators and function names:

\begin{hscode}
f :: Term Int -> Term Int
f x = x * x + 1
\end{hscode}

The library provides backend-specific implementations for many standard functions including all of the functions from such type classes as @Num@, @Floating@, @Ord@-like and @Eq@-like, as well as the ability to compose them into more complex functions.


\subsection{\name{HOAS} vs. \name{de Bruijn} representation}

Higher-order abstract syntax presented above is very convenient for library users to write. However, this representation is inconvenient for the embedded language compilers and interpreters to work with. In particular, it makes is impossible to inspect bodies of lambda functions.

The conversion of \name{HOAS} to an alternative representation based on \name{de Bruijn indices} has been independently discovered in \cite{HOAS2GADT} and \cite{Unembedding}. The conversion code employed in \LiveFusion is due to my supervisor Manuel Charkravarty, the author of the former approach.

The new representation offers an easier way to analyse embedded terms and is the assumed representation of scalar functions in the rest of this paper.



\subsection{Backend specific function implementations}

One notable part of the @Term@ language is its @Con@ constructor. The only argument to the constructor is @Impl t@ offers a backend-dependent way of to define a function. @Impl@ data type gives freedom to the backend to choose the most suitable representation for a function.

For instance, in the \Haskell backend for \LiveFusion we are interested in generating \Haskell source code, for which \name{Template Haskell} \cite{TH} is a natural choice. Thus we choose the following definition of @Impl@:

\begin{hscode}
data Impl t = HsImpl {
                hs :: t,        -- Native Haskell function
                th :: Q TH.Exp  -- TemplateHaskell quasiquoted expression
              } 
\end{hscode}


The \name{Template Haskell} expression is really the core part of the function representation, but this is what would later on allow for an easy production of \Haskell source code at the code generation stage.

The representation also includes the vanilla \Haskell function. At present this is just for completeness. However in the future it may be possible to avoid compiling certain parts of AST at runtime and run statically scheduled combinators. In this case this would be the function to inline into the loop.

We will now see the (rather trivial) implementations for a couple of functions before continuing with our explanation of Impl.

\begin{hscode}
plusImpl :: Num a => Impl (a -> a -> a)
plusImpl = HsImpl { hs = (+); th = [| (+) |] } 

absImpl :: Num a => Impl (a -> a)
absImpl = HsImpl { hs = abs; th = [| abs |] } 
\end{hscode}

It is easy to see that \name{Template Haskell} expressions are trivially created from simple \Haskell functions using \name{quasi-quotation} extension \cite{QQ}.

When the backend is ready to generate \Haskell source for these functions, it will be a simple matter of using a \name{Template Haskell's} pretty printer. 

The following are the reasons for using \name{Template Haskell} expressions to represent the user specified functions:
\begin{itemize}
\item The backend does non-trivial code generation as discussed in Chapter \ref{ch:Code-Generation}, where the full power of \name{Template Haskell} is required. Functions defined this way become very easy to incorporate into the generated code
\item New developments in \name{Template Haskell}\footnote{Starting with GHC 7.8.} will allow its expressions to be typed. Having @Q (TExp t)@ as opposed to @Q Exp@ will provide an extra layer of safety to the function representation
\end{itemize}

\todo{Perhaps show how a function is passed to the AST using lam and app.}
\todo{subsection{Function composition}}


%\subsection{User-defined implementation of parametrising functions}

%However, the flexibility of the library does not end there. One feature of the given scalar functions representation is that it allows the user of the library to define new functions. After all, it only requires creating a new record of type @Impl t@ where @t@ is the type of the desired function. Library user can then provide the internals of function implementation to appear in the generated code.

%This leaves the door open to the user and makes the library more flexible and extensible.


\clearpage
\section{Related work}

\subsection{\name{DESOLA}}
\label{sec:DESOLA}

\name{DESOLA} (\name{Delayed Evaluation Self Optimising Linear Algebra}) \cite{RMKB06} is an experimental \name{C++} library designed to explore the benefits of runtime code generation and optimisation for scientific computing. In contrast to \DPH, \name{DESOLA} does not support higher order combinators or segmented operations, although it does support regular multidimensional arrays.


\subsection{\name{Accelerate}}

\name{Accelerate} EDSL\iedsl \cite{CKL+11, McDonell:2013wi, CliftonEverest:2014vi} is probably the closest, by design, to \LiveFusion. It was initially designed to offer a simple but powerful programming model targetting General Purpose GPU computing. When I was starting my work on \LiveFusion, \name{Accelerate} was able to generate efficient code for GPUs but lacked fusion. Additionally, it only supported regular arrays (although with a very powerful \*shape-polymorphic* model similar to that of \Repa \cite{KCL+10}). \name{Accelerate} is currently moving towards having more sophisticated fusion, nested parallelism and ability to generate code for CPUs.


\subsection{\FlowFusion}
\label{sec:Fronted-FlowFusion}

Embedding EDSLs in a way that they construct their intermediate data structures at the runtime of the program is by far the most common approach but not the only one. An EDSL can also be implemented in a host language compiler.

In particular \GHC allows to specify transformations of its internal \name{Core} language in the form of plugins integrating into the compilation pipeline. \FlowFusion framework \cite{FlowFusion, FusingFiltersILP} was designed to fuse combinator graphs at compile time. The data flow graph to be fused is translated to \name{DDC Core} language -- the intermediate language of \name{Disciplined Disciple Compiler} \cite{DDC} -- where it is fused using  techniques that inspired the work on \LiveFusion and are discussed in the next chapter: rate inference (Section~\ref{sec:rates}) and compiling to a common imperative loop structure (Section~\ref{sec:anatomy}).

The approach of \FlowFusion solves sharing by analysing the entire graphs of combinators and avoids runtime overheads of building, analysis and compilation of ASTs. The downsides of \FlowFusion are that control flow breaks fusion. However, in \DPH control flow is not an issue since most of control flow is transformed into data flow by the vectoriser.

\FlowFusion is currently in its early development stages and is not yet able to fuse segmented combinators. In the future it may become a very comprehensive fusion platform. In fact there are plans for building the CPU backend for \Accelerate on top of \FlowFusion.



%% The following discussion is excluded from the document
\begin{comment}
\section{On the flexibility of parametrising functions}

\begin{itemize}
\item Currently function representation of only one backend was shown, namely Haskell source backend.
\item Even that one was trivial: only containing a TH Exp + native HS function
\item What other functionality may we want in our function representation?
\item LiveFusion is a fast array processing library. The focus is on performance.
\item Our two main options for speed on modern CPUs: parallelism + vector instructions.
\item Parallelism is already exploited by DPH through other means and will split the computations evenly across all cores
\item Vectorisation on the other hand has not been exploited by DPH yet\footnote{Though new developments in that area have emerged more recently \todo{ref}}
\item The CPU manufacturers are introducing and announcing more and more vector instruction sets\footnote{\todo{ref}: Intel Future Instruction Sets} which will allow many types of vector operations to be performed on each core.
\item This is very handy for libraries like LiveFusion which are already providing patterns of array computation, some of which map very well to vector instructions of current or future CPUs \todo{find examples, but this will definitely include point-wise arithmetic, reductions and pack}.
\item At some point in the future @Impl@ function representation may be the most natural place to include backend- or even CPU-specific information, on how to generate vector instructions for a particular function, e.g. a reduction using product operator on an Intel CPU.
\end{itemize}

\end{comment}

\begin{comment}
\section{Compiling Delayed LiveFusion AST}

\section{Sharing recovery: Abstract Semantic Graphs}

\begin{itemize}
\item Begin with Oleg's blessing quote. On two probs.
\item Find example, should be plenty in QuickHull
\item Talk about ref transparency of Haskell
\item Two types of sharing really, talk shapes
\item While can still fuse, it performs unnecessary work
\item Implicit vs. explicit sharing (ref email thread)
\item For implicit:
\item Still represented by the same object internally at runtime
\item StableNames, Gill's, but can do other, pointer equality sharing recovery
\item It's is IO, but since we are in a library and can guaranty safety by design we don't care
\end{itemize}
\end{comment}

\IfNotCompilingAll{\printbibliography}

\end{document}