\documentstyle{csjo}

\begin{document}
\bibliographystyle{ieeetr}

\title{A Grid-based Scientific Problem, Expressed in Haskell}

\author{D.B. Carpenter and H. Glaser}
\date{}
\maketitle

\section{Introduction}

     Lazy functional languages have not, to date, made an enormous
impact in scientific computing.  Partly, this has to do with
performance---codes written in these languages often run an order of
magnitude slower than imperative codes, which limits their appeal for
number-crunching production codes
(it seems that this caveat
no longer holds for some non-lazy functional
language---see \cite{scientific3}).
Apart from that, there is a
suspicion that programming without assignments or side-effects is
difficult, or restrictive.  Evidently, this is not the opinion of
functional programmers---but perhaps there are things in scientific
computation which really cannot be {\em done} easily without assignments.

     In a recent article \cite{self} we have given detailed Haskell
implementations of three ``scientific'' algorithms.
Haskell is a relatively new language, put forward as
an international standard for functional programming
\cite{haskell1,haskell2,haskell3}.  It is a lazy, pure functional language.
In order of increasing complexity, the problems considered were:
solution of a partial differential
equation by the relaxation method; Monte Carlo simulation of a
statistical mechanics system; and a similar simulation of a lattice
gauge theory.  This report describes a simplified version of
the first of those programs.
(Other examples of scientific functional programming can be found in
\cite{scientific1,scientific2,scientific3}).


\section{\label{poisson}The 2-dimensional Poisson equation}

     A fundamental problem in electrostatics is calculation of the
electric field produced by a given distribution of charges.  The
electrostatic potential is obtained by solving the associated Poisson
equation.  In two dimensions, the fixed
charge distribution is a function $c (x, y)$ of the coordinates.  The
electrostatic potential $u (x, y)$ is the solution of the equation
\begin{equation}
\label{continuum}
(\mbox{\boldmath $\nabla$}^2 u) (x, y) = - c (x, y).
\end{equation}

     The most simple-minded approach to solving (\ref{continuum}),
is to replace it by the finite difference approximation
\begin{equation}
\label{discrete}
\begin{array}[b]{l}
u (x + 1, y) + u (x - 1, y) + \\
u (x, y + 1) + u (x, y - 1) - 4 u (x, y) 
\end{array}
= c (x, y).
\end{equation}
The coordinates $x$, $y$ are now integers.  Rearranging (\ref{discrete}) as
\begin{equation}
\label{relaxation}
u (x, y) = \frac{1}{4} (\mbox{sum neighbouring $u$ values} + c (x, y)),
\end{equation}
immediately suggests the simplest iterative algorithm for solving this
system.  In the {\em relaxation} algorithm, starting from some initial guess
for $u$, each iteration replaces the value of $u$ at a point by the average
value of its neighbours (plus an offset determined by the local
charge).  Variants differ in whether the $u$ values inserted in the RHS
of (\ref{relaxation}) are those from the previous sweep of the grid (Jacobi
relaxation) or whether they are the most recent values, which may
include some values already updated in the current sweep
(Gauss-Siedel relaxation).  The last version has further variations
depending on the order in which the grid is traversed---a
popular scheme is ``red-black relaxation'' in which all even sites are
visited first, then all odd sites.


\section{Accessing neighbours}

     The two-dimensional geometry of the problem can be
embodied in the following definitions.
\begin{verbatim}
  type Coord   = (Int, Int)

  neighbs :: Coord -> [Coord]
  neighbs (x, y) = [(x, y + 1), (x, y - 1),
                    (x + 1, y), (x - 1, y)]
\end{verbatim}
\verb$neighbs$ returns a list of sites neighbouring its argument.
The remainder of the code will be independent of the lattice
dimension or shape (in \cite{self} we introduced some more
general techniques for achieving the goal of ``geometry independence'').
We define a specialised \verb$map$-like function for an array, which
operates on a pair-wise, nearest-neighbour basis
\begin{verbatim}
  nns :: (Coord -> Coord -> b -> c) ->
         Array Coord b -> Array Coord [c]
  nns f a =
    array bds [i :=
      [f i j (a ! j) | j <- neighbs i]
              | i <- range bds]
    where bds = bounds a
\end{verbatim}
(In Haskell, \verb$array$ is the array constructor function---it
takes an ``association list'', in which index values are
paired with data by the infix constructor \verb$:=$, and
returns an array.  This can be subscripted by \verb$!$.  It has
the same association between indices and data elements.)

\verb$nns$ is a function which takes a mapping function and
an array, and returns an array in which each element is
a mapped list of the nearest neighbour elements from the
original array.  If $j_1, j_2, ...$ are the neighbours of
site $i$, and \verb$b = nns f a$, then
\newline
\newline
\verb$  b ! $$i$\verb$ ==$
\newline
\verb$      [f $$i$\verb$ $$j_1$\verb$ (a ! $$j_1$\verb$), f $$i$\verb$ $$j_2$\verb$ (a ! $$j_2$\verb$), ...]$
\newline
\newline
Another higher order function we will use is \verb$azipWith$.
\begin{verbatim}
  azipWith f a b =
    array bds [i := f (a ! i) (b ! i) |
                i <- range bds]
    where bds = bounds a
\end{verbatim}
This is analogous to the standard \verb$zipWith$ function on lists.


\section{\label{implementations}Lazy implementation of Red-black
relaxation}

     With the geometry sorted out, we can tackle the main program.
Define a data type which captures the local
state in a relaxation calculation:
\begin{verbatim}
  type Obs    = Float

  data State  = Interior Float Obs |
                Boundary Obs

  siteObs :: State -> Obs
  siteObs (Interior _  u) = u
  siteObs (Boundary    u) = u
   
  siteTrans :: State -> [Obs] -> State
  siteTrans (Interior c _) nbs =
             Interior c ((sum nbs + c) /
                     fromIntegral (length nbs))
  siteTrans  b             _   = b
\end{verbatim}
There are two types of site: normal sites on the ``interior''
of the lattice, and those on some designated boundary (typically
the sites on the edge of the array).
Sites on the interior
have a local charge density and a potential value.  Those on the
boundary just have a potential value, which will remain fixed
throughout the iteration.  The function \verb$siteObs$
extracts the ``observable'', i.e. the potential, from a
\verb$State$ variable.  \verb$siteTrans$ performs the update of
equation (\ref{relaxation}) on the interior sites and leaves boundary
sites unchanged.  It is passed the old value of the site and a list
of potential values at nearest neighbours.

     The history of a single site can be regarded as a stream of
\verb$State$ values (by definition, a ``stream'' is
an infinite list).  The site evolves through application of
a transition function which depends on the (changing) environment of
the site.  This pattern of evolution can be captured through
the function.
\begin{verbatim}
  process :: (s -> e -> s) -> s -> [e] -> [s]
  process trans st ~(env : envs) =
      st : process trans (trans st env) envs
\end{verbatim}
(Note the lazy pattern-match on
the last parameter of {\bf process}.
This laziness is required for the implementations below.)
When the transition function is \verb$siteTrans$, the ``current environment''
is a list of neighbour
potential values, and the last argument to \verb$process$ will
be a stream of such lists---the history of the environment.

     The result of the computation will be an array of lists
of \verb$State$ values, defining the history of the whole lattice.
Suppose we are given this array and need to calculate
the array of environments from it.  This can be done by
\begin{verbatim}
  env :: Array Coord [State] ->
           Array Coord [[Obs]]
  env = amap transpose .
          nns (\_ _ -> map siteObs)
\end{verbatim}
The application of \verb$nns$, in which the mapping function
extracts a stream of potential values from the stream of
\verb$State$ values, gives an array containing a list of
streams of neighbour histories.  Mapping \verb$transpose$
over it turns these into streams of neighbour lists.
These are just what \verb$siteTrans$ needs as its argument.
So the implementation of relaxation is
\begin{verbatim}
  evolve :: Array Coord State ->
              Array Coord [State]
  evolve init = x
    where x = azipWith (process siteTrans)
                init e
          e = env x
\end{verbatim}
The argument \verb$init$ is an array of initial states (\verb$azipWith$ is
defined to pick up the lattice from its first array parameter).

     As well as the lattice size, the array of initial states defines:
which array elements are on the boundary of the physical system being
described; the charge distribution; the initial guess for the potential
on \verb$Interior$ sites; and the fixed values of the potential at
\verb$Boundary$ sites.  It's required that no \verb$Interior$ site
should be on the edges of the array (only \verb$Boundary$ sites are
allowed there).  Because \verb$siteTrans$ doesn't evaluate its neighbours
argument for \verb$Boundary$ sites, laziness ensures automatically
that no neighbour access violates array bounds.

     The code above implements Jacobi relaxation.
More general schemes, like Gauss-Siedel,
can easily be implemented in the same framework.
Suppose a site $A$ has site $B$ as a neighbour.
The scheme above calculates the new value at $A$ in terms
of the previous generation value at $B$.  To calculate
the value of $A$ in terms of the {\em latest} value of $B$,
we should look one place forward in the history of $B$---in
other words the $B$ stream should be tailed before
feeding it into the environment of $A$.
To organise this, it is convenient to attribute a colour
to each site.  If the colour of $B$, say, is less than
the colour $A$, $A$ is updated in terms of the latest
value of $B$.  Otherwise it is updated in terms of the
previous generation value of $B$.  The function \verb$env$
is changed to
\begin{verbatim}
  env = amap transpose .
          nns (\i j -> map siteObs .
            if colour j < colour i
              then tail
              else id)
\end{verbatim}
Red-black Relaxation in 2 dimensions can be implemented by setting
\begin{verbatim}
  colour (x, y) = (x + y) `mod` 2
\end{verbatim}


\section{\label{discussion}Discussion}

     The implementation creates an
``array of streams''.  An alternative approach is to create a ``stream
of arrays'', in which each update is performed on the whole lattice to
create a new global state.  It only involves minor changes to the
implementation to create such a stream of arrays
directly.  Alternatively one can feed the results of the given
implementation to a pipeline of \verb$elems$, \verb$transpose$ and
\verb$map (listArray bds)$ to produce a list of arrays.  There are some
advantages to the ``array of streams'' approach (apart from the fact
that it models process-parallelism).  Some existing Haskell
implementations do not handle array accesses very efficiently,
and in this
approach the array operations are confined to an initial stage where
the streams are ``tied together''.  In the ``stream of arrays''
approach array operations are needed in all updates.  On the other hand
it does seem more natural to have a list of arrays when it comes to
processing the data (the set of configurations) generated by a
simulation.

     The algorithm described involves state transition in an
essential way, but it is doubtful whether its description
would have been any clearer in terms of variables and assignments.
The explicit representation of sequences as infinite lists
seems just as natural.


\bibliography{haskell}

\end{document}
