
\documentstyle[harvard]{jfp}

\begin{document}
\bibliographystyle{agsm}

\sloppy

\title{Some Lattice-based Scientific Problems, Expressed in Haskell}

\author{D.B. Carpenter, H. Glaser, \\
        Department of Electronics and Computer Science, \\
        Southampton University, \\
        Southampton, S09 5NH} 

%\author{D.B. Carpenter}
%\address{Department of Electronics and Computer Science, \\
%         Southampton University, \\
%         Southampton, S09 5NH} 

%\date{DRAFT}

\maketitle

%\newpage

\begin{abstract}
The paper explores the application of a lazy functional
language, Haskell, to a series of grid-based scientific
problems---solution of the Poisson equation and Monte Carlo simulation
of two theoretical models from statistical and particle physics.  The
implementations introduce certain abstractions of grid topology,
making extensive use of the polymorphic features of Haskell.
Updating is expressed naturally through use of infinite lists,
exploiting the laziness of the language.  Evolution of systems
is represented by arrays of interacting streams.
\end{abstract}

%\tableofcontents
%\listoftables

%\newpage

%\setcounter{page}{1}
\pagestyle{plain}

%\input{introduction}

\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 at least an order of
magnitude slower than imperative codes, which limits their appeal for
number-crunching production codes.
(It seems that this caveat on
functional programming no-longer holds for some non-lazy functional
languages---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.

     This article gives detailed Haskell implementations of
three ``scientific'' algorithms.  In order of increasing
complexity they are: 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.  These problems are taken from theoretical
physics---a strongly mathematical discipline.  Functional
programming also has a very mathematical feel to it.
So it should not be too surprising that these problems actually
come out rather neatly in a functional language.
(Other examples of scientific functional programming can be found in
\cite{scientific1,scientific2,scientific3,scientific4,scientific5,scientific6}).

     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.  For the most part, the programs given here could have been
expressed just as well in any other lazy functional language.  Having
chosen Haskell, however, free use has been made of its more innovative
features, including the array and class systems (the latter for
achieving overloading, or {\em ad hoc} polymorphism, in the array
manipulation functions we introduce).  A fair understanding of Haskell
is probably a prerequisite for reading this article.  By and large, the
``prelude'' functions will be used without much comment, so a language
reference may be necessary.

     Full understanding of the problems and algorithms described will
require a certain amount of mathematical knowledge---some
familiarity with simple partial differential equations,
complex numbers, and a rough grasp of probability distributions.
Section \ref{problems} describes the systems involved
and the associated algorithms.  The problems are all based on lattices
(i.e., grids).  They make important use of arrays.

     Rather than work directly with standard Haskell arrays, it proved
convenient to introduce a kind of generalised array, with more
geometric significance.  Section \ref{lattices} describes the approach
taken to lattice geometry, and the arrays defined on these lattices.
The appendices give details of a possible Haskell implementation of the
generalised arrays and operations on them.

     The implementations of the main algorithms are given in section
\ref{implementations}.  The paradigm used in each case is one of
interacting parallel processes, modelled by lazy lists (streams)
representing the output of the processes.  Partly this choice can be
attributed to the source of funding for this work---the FAST project
aims to provide a parallel implementation of annotated functional
programs written in essentially this style \cite{fast1}.  Apart from
that, the style comes naturally and is quite efficient in a lazy
language, especially for algorithms with some ``locality of
reference''.

     Some conclusions are collected in section \ref{discussion}.

     We chose to handle the lattices in a rather abstract way.  The
result is that the geometry factors out from the problems in an
attractive way---the implementations in section \ref{implementations}
make almost no reference to the details of this geometry.  Abstraction
certainly has dividends, but understanding it also requires some
investment of effort.  Needless to say, Haskell does not {\em enforce}
this level of abstraction.  It does enable it.

%\input{problems}

\section{\label{problems}Three problems}

     The choice of problems here reflects an early motivation for this
work, which was to implement some QCD codes in a functional language.
{\em Quantum Chromodynamics} is the theory of the strong nuclear
force.  It has a reputation as a ``grand challenge'' numerical
problem.  In the present work (as it turned out) the emphasis is more
on elegant specification of the underlying algorithms than high
performance.

     Through a series of mathematical tricks
\cite{feynman,wilson1,wilson2,wilson3} QCD can be expressed as a {\em
lattice gauge theory}, in which form it is amenable to Monte Carlo
simulation.  In this article, lattice gauge theory will be introduced
from a more algorithmic point of view, assuming no particular knowledge
of the physical systems involved.  First two simpler problems which
share some of its features are introduced.

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

     A fundamental problem in electrostatics is calculating the
electric field produced by a given distribution of charges.  The
electrostatic potential is obtained by solving the associated Poisson
equation.  For simplicity we describe the two-dimensional case---a
third dimension introduces nothing essentially new.  The given
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}
Having solved for $u$, the electric field is the gradient
$\mbox{\boldmath $F$} (x, y) = - (\mbox{\boldmath $\nabla$} u) (x, y)$.

     The most simple-minded approach to solving (\ref{continuum}),
is to replace it by the finite difference approximation
\begin{equation}
\label{discrete}
u (x + 1, y) + u (x - 1, y) +
u (x, y + 1) + u (x, y - 1) - 4 u (x, y) = - c (x, y).
\end{equation}
The coordinates $x$, $y$ are now integers---assume that they and $c$
have been rescaled to eliminate appearance of the grid spacing parameter.
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-Seidel 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.  Straight relaxation is
a slow algorithm, but it does form the core update of more powerful
methods---see for example \cite{multigrid}.

\subsection{\label{xy}The XY Model}

     Statistical physics is largely concerned with the study of
phase transitions in complex systems.  Many simplified mathematical
models have been introduced which exhibit transitions similar to those
of physical systems.  The {\em XY model} is an example \cite{xyModel}.
This kind of model can be motivated by referring
back to the Poisson equation.  Equations (\ref{continuum}) or
(\ref{discrete}) can be obtained from variational principles.
Solving the system (\ref{discrete}) is equivalent to minimising the
functional
%\begin{equation}
%\label{electroEnergy}
%  E (u) = \sum_{x, y} \frac{1}{2} \left((u (x + 1, y) - u (x, y))^2 + 
%				        (u (x, y + 1) - u (x, y))^2\right) -
%	  c (x, y) u (x, y).
%\end{equation}
\begin{eqnarray}
  E (u) & = & \sum_{x, y} \frac{1}{2} \left((u (x + 1, y) - u (x, y))^2 + 
		      (u (x, y + 1) - u (x, y))^2\right) - \nonumber \\
\label{electroEnergy}
	&   & \hspace{1in} c (x, y) u (x, y).
\end{eqnarray}
$E (u)$ is the electrostatic energy of the configuration $u$; the
lowest energy state is the one ``preferred by nature''.

     In statistical physics, states with higher energy are also taken into
account.  The aim of a simulation is to generate a sequence of
configurations of the system, in which the frequency of appearance of a
particular configuration $u$ with energy $E (u)$ is weighted by the
Boltzmann probability distribution:
\begin{displaymath}
  P (u) \propto \exp (- E (u) / \theta).
\end{displaymath}
Here, $\theta$ is the temperature of the system.

     In the XY model $u (x, y)$ is a unit-length two-dimensional vector
rather than a simple real number.  This vector can be represented as a
complex number with unit magnitude.  The electrostatic energy function,
(\ref{electroEnergy}), is replaced by
\begin{equation}
\label{XYEnergy}
  E (u) = - \Re \sum_{x, y} u^* (x, y) u (x + 1, y) +
                            u^* (x, y) u (x, y + 1),
\end{equation}
where $u^*$ means complex conjugate of $u$, and $\Re$ means real part.

     Algorithms for generating sequences of configurations with Boltzmann
distribution in this type of system resemble relaxation.  They are
iterative algorithms, involving repeated sweeps across the grid, and the
update of a single site normally only needs knowledge of the current
values at neighbouring sites.

     In the {\em Metropolis} algorithm \cite{metropolis} 
the update of a site is as
follows:  choose a candidate new value for the local $u$ by some
pseudo-random procedure.  The procedure used for this choice should give equal
probabilities for rotating the vector by $+\delta$ and $-\delta$, for
all angles $\delta$.  Also it should be ``ergodic''---by combining
some sequence of allowed updates it must be possible to get
from any state to any other.  Calculate the change in energy, $\Delta E$,
which would be produced by replacing the old value of $u$ by the
candidate new one (this depends on the state of the neighbouring
sites).  If $\Delta E$ is negative, accept the candidate update.  If
$\Delta E$ is positive, accept the candidate with probability $P = exp
(-\Delta E / \theta)$.  An operational procedure for doing this is to
generate a pseudo-random value $x$ uniformly in the interval $[0, 1)$,
accepting the candidate if $x < P$.  If the candidate is rejected,
retain the original value of $u$.

The proof of this algorithm is omitted.  The crux
is to show that, when the probability of a system state is given by the
Boltzmann formula, the expected frequency of transitions
(produced by the algorithm) from any state $A$ to any other state $B$
is the same as that from $B$ to $A$ (``detailed balance'') so
that the Boltzmann distribution represents equilibrium.

\subsection{\label{lgt}Lattice gauge theory}

     This section describes the so-called $U(1)$ lattice gauge
theory.  $U(1)$ is the group of unit magnitude complex numbers,
with complex multiplication as the product.  This is the ``gauge
group'' of {\em Quantum Electrodynamics} rather than QCD, but
the two theories are closely related.
\begin{figure}[hbt]
\centerline{
\begin{picture}(80, 90)(-20, -30)
\multiput(-10,  0)( 0, 20){3}{\line(1, 0){60}}
\multiput( 0, -10)(20,  0){3}{\line(0, 1){60}}
\put(0, 0){\circle*{6}}
\put(10, -30){Site}
\end{picture}
\begin{picture}(80, 90)(-20, -30)
\multiput(-10,  0)( 0, 20){3}{\line(1, 0){60}}
\multiput( 0, -10)(20,  0){3}{\line(0, 1){60}}
\put(2,  0){\rule[-1pt]{16pt}{3pt}}
\put(10, -30){Link}
\end{picture}
\begin{picture}(80, 90)(-20, -30)
\multiput(-10,  0)( 0, 20){3}{\line(1, 0){60}}
\multiput( 0, -10)(20,  0){3}{\line(0, 1){60}}
\put(2,  2){\rule{16pt}{16pt}}
\put(0, -30){Plaquette}
\end{picture}
}
\caption{\label{cells}Some elementary cells of a lattice.}
\end{figure}
     The main difference between $U(1)$ gauge theory and the XY
model is geometric.  In the XY model the variables $u$ are
associated with the {\em sites} of the lattice.  The energy
function (\ref{XYEnergy}) is a sum over the {\em links} of the
lattice---it involves multiplying together the values of the
variables at the ends of each link.

     In gauge theory the variables are associated with the links
of the lattice, and the energy function is computed on the
next-higher dimension cells of the lattice---the elementary
squares, or {\em plaquettes} (figure \ref{cells}).

     We can associate a group element (a unit-magnitude complex number)
$U_P$ with any path $P$ in the lattice by multiplying 
together the link variables along the path.  In particular, we
can associate a path with each elementary plaquette---the four-link
path bounding that plaquette.  The energy function is now
\begin{displaymath}
E (u) = - \sum_{P \in \mbox{\em plaquettes}} \Re (U_P)
\end{displaymath}
where {\em plaquettes} is the set of elementary plaquette-paths in the
lattice.  (In the XY model energy, one of the factors in the link product
has to by conjugated.  In gauge theory, variables associated
with links that have negative incidence on the plaquette are conjugated.
This will be illustrated more fully in the next two sections.)

     The Metropolis algorithm can be carried over from the XY model.
The local update is identical.  The difference lies in how the
neighbouring values that affect the local contribution to the
energy are collected.  In the XY model, the local contribution
to the energy was determined by the values of site variables
at neighbouring sites.  In gauge theory it is determined
by link variables on the ``staples'' illustrated (in the 3-dimensional
case) in figure \ref{staples}.  For real physics, the lattice should
be four-dimensional.  It is usually taken as a 4d cubic grid, but the
theory has sometimes been formulated on other lattice geometries;
the implementations given here will be largely independent of
the details of this geometry.
\begin{figure}[hbt]
\centerline{
\begin{picture}(120, 100)(0, -20)
\thicklines
\put(60, 0){\line(0, 1){60}}
\thinlines
\put(65, 5){\line(1, 1){20}}
\put(65, 65){\line(1, 1){20}}
\put(85, 25){\line(0, 1){60}}
\put(55, -5){\line(-1, -1){20}}
\put(55, 55){\line(-1, -1){20}}
\put(35, -25){\line(0, 1){60}}
\put(0, 0){\line(0, 1){60}}
\put(0, 0){\line(1, 0){50}}
\put(0, 60){\line(1, 0){50}}
\put(120, 0){\line(0, 1){60}}
\put(70, 0){\line(1, 0){50}}
\put(70, 60){\line(1, 0){50}}
\end{picture}
}
\caption{\label{staples}Surrounding links, contributing to the energy of the central one.}
\end{figure}

%\input{lattice}

\section{\label{lattices}Generalised Arrays}

     A straightforward way to store the values of, for example,
the potential function $u$ in section \ref{poisson} is in
a Haskell array.  This works well, provided the lattice is a
regular grid and its extent is rectangular.  It works less
well in the case of lattice gauge theory (section \ref{lgt}).
There the variables naturally live on the links of the lattice.
One can still use an array indexed by lattice sites, with
an extra dimension taking values in $1, \ldots, d$ to access the
$d$ links rooted at each site ($d$ is the dimension of the lattice).
This is quite clumsy.

     We can think of the sites of a lattice as the 0-dimensional
cells of the lattice, the links as the 1-dimensional cells, the
plaquettes as 2-dimensional cells, and so on.  Occurrence
of variables which naturally live on cells with
dimension higher than 0 is more common than one might initially
suppose.  Even in the Poisson equation the electric field, which
is the gradient of the potential, sits most naturally on the
links of the lattice.  The field in the $x$-direction
is associated with links in that direction; its value is just the
difference of the potential at the two ends of its link.
In the XY model there occur certain ``vortex configurations'' centred
on the plaquettes of the lattice---one can define a function which
counts the density of vortices in a general configuration, and this
function sits naturally on the plaquettes of the lattice.  As
explained in section \ref{lgt}, the energy function of a gauge theory
also ``lives'' on the plaquettes of a lattice.

     This suggests that it may be useful to generalise the usual model
of an array so that it is defined on an arbitrary (not necessarily
rectangular) lattice, and that the index space of a particular array can
be cells of that lattice of arbitrary (fixed) dimension.  This
generalised kind of array will be called a ``form''.  The
terminology is motivated by an analogy between this kind of array and
the $p$-forms of differential geometry \cite{geometry}.

     The lattice itself is no longer thought of just as a set of
connected sites, but as a ``cell complex'' containing cells of
all dimensionality up to that of the lattice itself.

\subsection{Cells, lattices and forms}

     This section assumes some knowledge of the Haskell class
system.  The standard prelude defines, for example, the class
\verb$Ix$.  All types which can be used as index types for an
array belong to this class.  A type specification involving
type variables may be preceded by an assertion that some of
those variables belong to class \verb$Ix$.
Array operations like the array subscripting operator \verb$(!)$
(which, by the way, is an ordinary infix function---not special
syntax) are {\em overloaded} to apply to all instances of
\verb$Ix$.

     We need a Haskell class \verb$Cell$ which parallels the
index class \verb$Ix$ associated with ordinary arrays.
Two type constructors \verb$Lattice$ and \verb$Form$ go with this class.
They construct the following types:
\begin{verbatim}
  (Cell a) => Lattice a
  (Cell a) => Form a b
\end{verbatim}
The \verb$(Cell a) =>$ prefix asserts that the type \verb$a$
must be an instance of class \verb$Cell$.
The second type above is directly analogous to the ordinary Haskell
\verb$Array$ constructor---\verb$Form a b$ is the type of a form with
cell-type \verb$a$ and elements of type \verb$b$.  \verb$Lattice a$
is the type of a data structure describing a lattice with
cells of type \verb$a$.  It takes over the role of the
bounds pair for an ordinary array.

     To make this more concrete, a possible instance of the
\verb$Cell$ class, suitable for handling 2-dimensional rectangular
lattices, is
\begin{verbatim}
  data Cell2d = Cell2d0 (Int, Int)     |
                Cell2d1 (Int, Int) Int |
                Cell2d2 (Int, Int)
\end{verbatim}
which can represent a site, link or plaquette of the lattice.
\verb$Form$ and \verb$Lattice$ types will be treated, in this section
at least, as abstract types.  The cell types could also be made
abstract---with the usual advantage of hiding implementation and
leaving room for improving it without changing code that uses the data
types.  One would introduce true functions (lower case) for
constructing the cell objects, and a set of access functions for the
coordinates, etc.  This would probably be sensible practice---forgone
here to save space on definitions.

The numeric fields are the
cell coordinates, and in the 1-cell case, the link direction.
The cell type alone does not
define the lattice shape.  We need at least to know the size of a
rectangular lattice, and besides that, the ``edge conditions'' have to
be specified---the lattice may have a boundary at the edges, or be
connected cyclically, etc.  We will assume the existence of
a predefined function such as
\begin{verbatim}
  torus2d :: Int -> Int -> Lattice Cell2d
\end{verbatim}
which takes the width and height of the lattice as parameters, and
returns an instance of an abstract data type representing (in this
example, let's say) a periodic lattice of the specified size.
A possible implementation of \verb$torus2d$ is given in appendix
\ref{topologyExample}.  Now, we can
start building forms.  The function analogous to the ordinary
Haskell \verb$array$ is
\begin{verbatim}
  form :: (Cell a) => Lattice a -> Int -> [Assoc a b] -> Form a b
\end{verbatim}
This constructs a form from an ``association list''.
The first parameter is the lattice on which the form is to be
defined; the second is the ``rank'', $r$, of the form---that is, the
dimension of the cells on which it is defined; the third is
an association list completely analogous to the corresponding
argument of \verb$array$---a list of pairs associating each cell
of dimension $r$ in the lattice with a value.  As a mundane
example, the expression
\begin{verbatim}
  form l 1 [c := 0.0 | c <- cells l 1]
\end{verbatim}
is 1-form with value $0.0$ on each link of lattice \verb$l$.
\begin{verbatim}
  cells :: (Cell a) => Lattice a -> Int -> [a]
\end{verbatim}
is analogous to \verb$range$, returning a list
of all cells of a specified dimension in the given lattice.
Ultimately, the following operations on forms will also be needed
\begin{verbatim}
  accumForm :: (Cell a) =>
    Lattice a -> Int -> [Assoc a b] -> Form a [b]
  lattice   :: (Cell a) => Form a b -> Lattice a
  rank      :: (Cell a) => Form a b -> Int
  (!#)      :: (Cell a) => Form a b -> a -> b
 
  (//#)     :: (Cell a) =>
    Form a b -> [Assoc a b] -> Form a b
  
  listForm  :: (Cell a) =>
    (Lattice a) -> Int -> [b] -> Form a b
 
  fassocs   :: (Cell a) => Form a b -> [Assoc a b]
  findices  :: (Cell a) => Form a b -> [a]
  felems    :: (Cell a) => Form a b -> [b]
 
  fmap      :: (Cell a) =>
    (b -> c) -> Form a b -> Form a c
  fmapi     :: (Cell a) =>
    (a -> b -> c) -> Form a b -> Form a c
  fzipWith  :: (Cell a) =>
    (b -> c -> d) -> Form a b -> Form a c -> Form a d
\end{verbatim}
Most of these are fairly self-explanatory generalisations of similarly-named
Haskell prelude functions.  A possible implementation of them
is given in appendix \ref{formLibrary}.  Use of those
definitions requires an instance of the cell class.  An example
instance is given in appendix \ref{cellInstance}.

     \verb$(!#)$ is the infix subscripting operation on a form, so
if \verb$a$ is a form and \verb$c$ is a cell of the right type and
dimension, \verb$a !# c$ is the corresponding element of \verb$a$.
\verb$lattice$ and \verb$rank$ generalise the standard \verb$bounds$
function, returning the lattice and rank of a given form.
\verb$accumForm$ is similar to \verb$form$ but multiple
associations are allowed for individual cells---the result is
a form containing, for each cell, a list of all values the association
list binds to that cell.
\verb$(//#)$ is for partial updates---it also resembles \verb$form$ but
takes a form as input.  The result is a new form with the same lattice
and rank; any cells unspecified in the association list take their
values from the input form.
\verb$fmap$ is a straightforward
mapping function.  \verb$fmapi$ is similar, but the mapping function
is passed an extra argument, the cell on which
the function map is being applied.  \verb$fzipWith$ is similar
to the \verb$zipWith$ function on lists.  It takes two conforming
forms and produces a third by applying the zipping function
at each cell.

\subsection{\label{latticeTopology}Lattice topology}

     The final and most interesting ingredient is the \verb$boundary$
operation.

     A powerful way of representing the topology of a cell complex
is through its {\em incidence function}.  This is a function $I$
from pairs of cells to the set $\{-1, 0, +1\}$.  Its value is
non-zero only when the cell pair is of the form $(c^{(r)}, c^{(r - 1)})$,
where $c^{(r)}$ is an $r$-dimensional cell, and $c^{(r - 1)}$
is an $(r-1)$-dimensional cell, {\em and $c^{(r - 1)}$ is on the
boundary of $c^{(r)}$}.  When it is non-zero, its value is $+1$ or
$-1$ according to whether $c^{(r - 1)}$ is positively or negatively
oriented on the boundary of $c^{(r)}$.  Orientation is easier to
grasp through pictures than by formal definition---see
figure \ref{incidence}.
\begin{figure}[hbt]
\centerline{
\begin{picture}(220, 130)(-20, -30)
\multiput(-10,  0)( 0, 40){3}{\line(1, 0){100}}
\multiput( 0, -10)(40,  0){3}{\line(0, 1){100}}
%
\put(40, 80){\circle*{6}}
\put(30, 84){$H$}
\put(80, 80){\circle*{6}}
\put(84, 70){$G$}
\put(45, 80){\rule[-1pt]{30pt}{3pt}}
\put(57, 70){$F$}
%
%\put(5,  5){\rule{30pt}{30pt}}
\put(5,  5){\rule{10pt}{30pt}}
\put(5,  5){\rule{30pt}{10pt}}
\put(5, 26){\rule{30pt}{10pt}}
\put(26,  5){\rule{10pt}{30pt}}
\put(17, 17){$A$}
%
\put(5,  0){\rule[-1pt]{30pt}{3pt}}
\put(17, -10){$B$}
\put(5, 40){\rule[-1pt]{30pt}{3pt}}
\put(17, 44){$D$}
\put(-1,  5){\rule{3pt}{30pt}}
\put(-10, 17){$E$}
\put(39,  5){\rule{3pt}{30pt}}
\put(44, 17){$C$}
%
\put(100, 110){
\begin{minipage}[t]{2in}
\begin{eqnarray*}
I (A, B) & = & +1 \\
I (A, C) & = & +1 \\
I (A, D) & = & -1 \\
I (A, E) & = & -1 \\
I (F, H) & = & -1 \\
I (F, G) & = & +1 \\
I (A, F) & = & 0 \\
I (B, H) & = & 0
\end{eqnarray*}
\end{minipage}
}
\end{picture}
}
\caption{\label{incidence}Example values of the incidence function on a 2d lattice.}
\end{figure}
In the figure the conventions are adopted that the links are directed
from left to right or from bottom to top, and the boundary of
a plaquette is oriented anti-clockwise around the square.
So, for example, link $B$ has positive incidence on plaquette $A$,
because it follows the direction of the boundary path, whereas $D$ has
negative incidence because it is directed left to right while
the top boundary of $A$ goes right to left.  For the incidence
of a site on a link, the convention is that it is positive
if the link is directed into the site, and negative if it
is directed out of the site.  Reversing any of these conventions on
orientation may change local signs in the incidence function, but {\em not} the
final, scalar results of computations.

     Now the \verb$boundary$ function,
\begin{verbatim}
  boundary  :: (Cell a, Num b) => Lattice a -> a -> [Assoc a b],
\end{verbatim}
applied to an $r$-cell \verb$c$
returns the list of $(r-1)$-cells on the boundary of that cell,
paired with the incidence function of the boundary cell on \verb$c$.
Taking examples in figure \ref{incidence}:
\begin{tabbing}
\verb$  $\= \kill
\>\verb$boundary l $$A$\verb$ = [$$B$\verb$ := 1, $$C$\verb$ := 1, $$D$\verb$ := -1, $$E$\verb$ := -1]$ \\
\>\verb$boundary l $$F$\verb$ = [$$H$\verb$ := -1, $$G$\verb$ := 1]$ \\
\>\verb$boundary l $$H$\verb$ = []$
\end{tabbing}
where \verb$l$ is the associated lattice data structure.
The result for $H$ is empty, because by convention 0-cells
have no boundary.
(There is no huge significance in the choice of \verb$Assoc$
pairs as the result of \verb$boundary$.  In principle, any type capable of
representing a cell and a numeric object would do.)
Appendix \ref{topologyExample} shows how a particular lattice
topology can be set up---defining an instance of the
\verb$boundary$ function.

     The \verb$boundary$ operation probably looks a bit contrived
at first sight.  The idea comes from homology---see for example \cite{topology}.
As the examples in the next section show, it is a convenient
practical method of encoding the shape of a lattice.

    To close this section we record a couple more functions which will be
used later, and which occur frequently in other applications of the
``form'' technology introduced above.  First, a function which
``contracts'' a form on an association list containing cells and values
of the same type:
\begin{verbatim}
  cntrct :: (Cell a, Num b) => Form a b -> [Assoc a b] -> b
  cntrct a c = sum [v * (a !# j) | j := v <- c]
\end{verbatim}
(This is a sort of scalar product between a form and an association
list.)

    Second, \verb$mltply$ multiplies all values
in an association list by its first argument.
\begin{verbatim}
  mltply :: (Cell a, Num b) => b -> [Assoc a b] -> [Assoc a b]
  mltply m ivs = [i := m * v | i := v <- ivs]
\end{verbatim}

%\input{implementations}

\section{\label{implementations}Form-based implementation of the three algorithms}

\subsection{Accessing neighbours}

Consider the function
\begin{verbatim}
  neighbours :: (Cell a, Num b) =>
                Lattice a -> Int -> Form a [[Assoc a b]]
  neighbours l r = accumForm l r
    [j := mltply v b | b <- [boundary l i | i <- cells l (r + 1)],
                       j := v <- b]
\end{verbatim}
\verb$neighbours l r$ is an $r$-form on the lattice $l$.  We claim that
the element of this form associated with a cell $j$ contains a certain
representation of the set of $r$-cells neighbouring $j$.

The above claim will be illustrated in
the case where $r$ is zero.  The list-comprehension variable
\verb$b$ is instantiated in turn to the boundary of every link in the
lattice.  Link $B$ in figure \ref{neighbours} would yield (up to a
possible re-ordering of lists)
\begin{tabbing}
\verb$  $\= \kill
\>\verb$b = [$$A$\verb$ := -1, $$F$\verb$ :=  1]$
\end{tabbing}
\begin{figure}[hbt]
\centerline{
\begin{picture}(220, 130)(-20, -30)
\multiput(-10,  0)( 0, 40){3}{\line(1, 0){100}}
\multiput( 0, -10)(40,  0){3}{\line(0, 1){100}}
%
\put(40, 40){\circle*{6}}
\put(30, 30){$A$}
\put(45, 40){\rule[-1pt]{30pt}{3pt}}
\put(57, 44){$B$}
\put(39, 45){\rule{3pt}{30pt}}
\put(44, 57){$C$}
\put( 5, 40){\rule[-1pt]{30pt}{3pt}}
\put(17, 44){$D$}
\put(39,  5){\rule{3pt}{30pt}}
\put(44, 17){$E$}
\put(80, 40){\circle*{6}}
\put(70, 30){$F$}
\put(40, 80){\circle*{6}}
\put(30, 70){$G$}
\put( 0, 40){\circle*{6}}
\put(-10,30){$H$}
\put(40,  0){\circle*{6}}
\put(30,-10){$I$}
%
\put(100, 110){
\begin{minipage}[t]{2in}
\begin{eqnarray*}
I (B, A) & = & -1 \\
I (C, A) & = & -1 \\
I (D, A) & = & +1 \\
I (E, A) & = & +1 \\
I (B, F) & = & +1 \\
I (C, G) & = & +1 \\
I (D, H) & = & -1 \\
I (E, I) & = & -1
\end{eqnarray*}
\end{minipage}
}
\end{picture}
}
\caption{\label{neighbours}Collection of 0-cells neighbouring $A$.}
\end{figure}
This gives rise to the terms
\begin{tabbing}
\verb$  $\= \kill
\>\verb$[ ..., $$A$\verb$ := [$$A$\verb$ :=  1, $$F$\verb$ := -1],$ \\
\>\verb$       $$F$\verb$ := [$$A$\verb$ := -1, $$F$\verb$ :=  1], ...]$
\end{tabbing}
in the outer list comprehension (the list in the second association
is obtained after multiplication by $-1$).
After the application of \verb$accumForm$, the value
at site $A$ will be (again, up to possible re-ordering of lists)
\begin{tabbing}
\verb$  $\= \kill
\>\verb$[[$$A$\verb$ :=  1, $$F$\verb$ := -1],$ \\
\>\verb$ [$$A$\verb$ :=  1, $$G$\verb$ := -1],$ \\
\>\verb$ [$$A$\verb$ :=  1, $$H$\verb$ := -1],$ \\
\>\verb$ [$$A$\verb$ :=  1, $$I$\verb$ := -1]]$
\end{tabbing}
The form produced by \verb$neighbours$ contains some redundant
information (probably all that was really needed was
the identities of $F$, $G$, $H$ and $I$).
However, the definition here is generalised
to produce forms of neighbours for arbitrary-dimension lattice cells.
This will be useful later.

     Now we define a mapping function for the 0-cell case, which
applies a map to neighbours of each cell
\begin{verbatim}
  nns :: (Cell a) => (a -> a -> b -> c) -> Form a b -> Form a [c]
  nns f a = fmapi (\i ivss ->
              [f i j (a !# j) | ivs <- ivss, j := v <- ivs, v < 0])
                  (neighbours l 0)
            where l = lattice a
\end{verbatim}
This exploits the knowledge that neighbours proper appear with a
negative value in the association list, filtering out occurrences of
the ``destination'' cell itself.  The map function is passed the
destination cell, the ``source'' cell (a neighbour to the destination
cell) and the value of the original form at the source cell.  In the
example above we would get, for ``destination'' $A$,
\begin{tabbing}
\verb$  $\= \kill
\>\verb$(nns f a) !# $$A$\verb$ = [$\=
    \verb$f $$A$\verb$ $$F$\verb$ (a !# $$F$\verb$),$
    \verb$f $$A$\verb$ $$G$\verb$ (a !# $$G$\verb$),$ \\
\>\>\verb$f $$A$\verb$ $$H$\verb$ (a !# $$H$\verb$),$
    \verb$f $$A$\verb$ $$I$\verb$ (a !# $$I$\verb$)]$
\end{tabbing}

\subsection{\label{relaxImplementation}Poisson equation relaxation.}

     With the geometry well sorted out, we can turn to our first
application.  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@(Boundary _) _   = b
\end{verbatim}
There are two types of site: normal sites on the ``interior'' of the
lattice, and those on some designated boundary (which could be empty).
This grid ``boundary'' should not be confused with the cell
\verb$boundary$ operation!  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 neighbouring sites.

     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
a function with the property
\begin{verbatim}
  process :: (s -> e -> s) -> s -> [e] -> [s]
  process trans st [e1, e2, ...] =
    [st, trans st e1, trans (trans st e1) e2, ...]
\end{verbatim}
\verb$st$ is the initial value of the state, \verb$e1$ is the
initial value of the environment, \verb$e2$ is the next value of
the environment, and so on.  This function already exists in
the Haskell prelude; we can set
\begin{verbatim}
  process = scanl
\end{verbatim}
In the case at hand 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 a 0-form of lists
of \verb$State$ values, defining the history of the whole lattice.
Suppose we are given this form and need to calculate
the 0-form of environments from it.  This can be done by
\begin{verbatim}
  env :: (Cell a) => Form a [State] -> Form a [[Obs]]
  env = fmap 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 a form 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 :: Form Cell2d State -> Form Cell2d [State]
  evolve init = x where x = fzipWith (process siteTrans) init e
                        e = env x
\end{verbatim}
The argument is a 0-form of initial states (\verb$fzipWith$ is
defined to pick up the lattice from its first array parameter).

     This code implements Jacobi relaxation.
More general schemes, like Gauss-Seidel,
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 = fmap 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
\begin{verbatim}
  colour (Cell2d0 (x, y)) = (x + y) `mod` 2
\end{verbatim}

\subsection{\label{XYimplementation}XY model simulation}

     To describe the local state in an XY model simulation,
the \verb$State$ data type is changed in the following way:
\begin{verbatim}
  type Obs   = Complex Float

  data State = MkStat Float [Float] Obs
 
  siteObs (MkStat _ _ s) = s
 
  siteTrans (MkStat theta rans spin) neighbs =
    MkStat theta rans'' spin'
    where
      (cand, rans')   = metroCand rans spin
      eDiff           = deltaE neighbs spin cand
      (spin', rans'') = metroChoice spin cand (eDiff / theta) rans'
\end{verbatim}
The fields in the \verb$State$ record are respectively the temperature,
an infinite list of random real numbers, and the local value of $u$ (sometimes
called the ``spin'').  Each site must hold an {\em independent}
stream of random numbers.  Creating the random number streams is part
of initialisation---we omit details here (see appendix \ref{running}),
but note that there are various approaches to getting independent
streams of pseudo-random numbers.  If the basic generator has a long
enough cycle, simply taking independent seed values, perhaps chosen by
a different generator, may be acceptable.  Alternatively, there are
various ``leapfrog'' schemes available---see, for example
\cite{leapfrog}---which are quite elegant.  Another attractive
approach, applicable even when the number of streams is not known a
priori, is described in \cite{burton}.  We actually used a leapfrogged
64-bit linear congruential generator.

     As explained in section \ref{xy}, the Metropolis update has three
stages---choose a candidate update; calculate the change it
would make to the system's energy; then decide on the basis of this
whether to accept the candidate or retain the old value at the
site.  The details are:
\begin{verbatim}
  -- Energy difference between two spin states.
  deltaE neighbs s s' =
      - realPart (conjugate (sum neighbs) * (s' - s))
 
  -- Candidate new spin for Metropolis update.
  metroCand (r : rans) old = (exp (0 :+ delta) * old, rans)
                             where delta = max_delta * (2 * r - 1)
 
  -- Metropolis choice.  Inputs are old and candidate new variable,
  -- energy difference of configurations, and a list of randoms.
  metroChoice _   new ediff  rans | ediff < 0 = (new, rans)
  metroChoice old new ediff (r : rans') =
      (if accept then new else old, rans')
      where accept = r < exp (- ediff)
\end{verbatim}

     \verb$metroCand$ returns a unit-magnitude 2-vector which differs
from the one input by rotation through a random angle in the range
$[ -$\verb$max_delta$$, +$\verb$max_delta$$ )$.  Here \verb$max_delta$
is a global variable---that would not be very satisfactory treatment in a
practical program, but it avoids some clutter here.

     The remainder of the code is identical to that in the previous
section.  The Red-black variant is preferred---in this case the
``Jacobi'' version would actually be invalid because neighbouring sites
connected by the energy function would be updated ``concurrently'', and
localised calculation of the energy change would be erroneous.

     In this and the preceding section we have shown how to generate an
array containing the histories of individual sites, without discussing
how the initial state is constructed, or how the output information is
typically used.  Appendix \ref{running} fills in some of these details.


\subsection{\label{lgtImplementation}Lattice gauge theory simulation}

     The 1-form analogue of \verb$nns$ is
\begin{verbatim}
  staples :: (Cell a, Num b) =>
    (a -> [(a, b, c)] -> d) -> Form a c -> Form a [d]
  staples f a =
    fmapi (\i ivss ->
           [f i [(j, v, a !# j)
               | j := v <- ivs, j /= i] | ivs <- ivss])
          (neighbours l r)
    where l           = lattice a
          r           = rank a
\end{verbatim}
The mapping function has a more complicated specification here.  There
are two reasons for this.  First, each adjacent plaquette
defines several (in a cubic lattice, three) neighbours to a given link,
whereas an adjacent link only defined a single neighbour to a given
site.  Second, orientation information is irrelevant to sites, but
important for links (links are directed; sites are not).  The mapping
function again takes the ``destination'' cell as its first argument.
Its second argument is a list.  The items in this list correspond to
the links in a single staple (see figure \ref{staples}).  Each
item is a tuple containing the coordinate (cell) of the link, its
orientation relative to the destination cell ($+1$ or $-1$ depending on
whether it points in the same or the opposite direction
around the boundary of the plaquette), and the
value held at the link.

     Now, the definition of \verb$env$ is changed to
\begin{verbatim}
  env = fmap transpose .
        staples (\i lss ->
	  (map product . transpose)
	  [((if v == -1 then map conjugate else id) .
	    map siteObs .
	    (if colour j < colour i then tail else id)) ls |
					  (j, v, ls) <- lss])
\end{verbatim}
The list comprehension pre-processes the streams of states associated
with the links on a given staple.  In contrast to the previous
version, variables are conjugated here if they have negative orientation
relative to the destination link (in the XY model, this conjugation was
also necessary, but was deferred to \verb$deltaE$).  The
resulting streams are zipped together by the
\verb$transpose$, \verb$map product$
pipe, to produce a stream of products around the staple.

     Colouring is slightly more complicated than before.
It is necessary to ensure that no two links lying on the same
plaquette have the same colour.  For a 4 dimensional cubic array
a suitable definition might be
\begin{verbatim}
  colour (Cell4d1 (t, x, y, z) d) =
      4 * ((t + x + y + z) `mod` 2) + d
\end{verbatim}

     To save space, this article does not define an implementation of
4-dimensional lattice geometry.  The modifications necessary to the
2-dimensional case given in the appendices are reasonably
straightforward.

     Finally, the local state and updating can be lifted verbatim from
the XY model, except that \verb$deltaE$ is simplified to
\begin{verbatim}
  deltaE neighbs s s' = - realPart (sum neighbs * (s' - s))
\end{verbatim}
because the conjugation was already done in \verb$env$.
The top-level evolution is unchanged from section \ref{relaxImplementation}
(the type specification of \verb$evolve$ has to be altered
slightly).  One might want to change some names: the ``site''
prefixes ought to become ``link'' or, more generically, ``cell''.

     Incidentally, we have been describing a theory with an Abelian gauge
group.  Minor changes are needed to \verb$neighbours$ and \verb$staples$
to support non-Abelian theories (QCD), because the order of the product
around the staples becomes important.  Hence, one has to take a little
more care to preserve ordering of the plaquette boundary lists.

%\input{discussion}

\section{\label{discussion}Discussion}

     The algorithms described involve state transition in an essential
way, but it is doubtful whether their description would have been any
clearer in terms of variables and assignments.  The explicit
representation of sequences as infinite lists seems just as natural.
In some ways the list representation is closer to traditional
mathematical notation, where one typically represents state
transitions by introducing a subscripted sequence of states.

     The implementations in section \ref{implementations} create 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$felems$, \verb$transpose$ and
\verb$map (listForm l r)$ 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.  Appendix \ref{running} ends with a simple example
of this kind of data processing.

     In any case, once the ideas about lattices and arrays in section
\ref{lattices} are absorbed, implementation of the three
algorithms is very straightforward.  It was worth the effort of doing
the geometry carefully, for the level of generality it gave in
the implementations.  For example, apart from the definition of
\verb$colour$, the modifications made to the XY model simulation in
section \ref{lgtImplementation} to encompass lattice gauge theory
did not actually invalidate the code as an implementation of the XY
simulation---the final version works for systems with variables defined
on arbitrary dimension cells of an arbitrary dimension lattice (given
an appropriate instance of the \verb$Cell$ class).
Achieving this level of abstraction depended on
both ordinary (parametric) polymorphism and overloading,
together with use of higher order functions.

     The price one pays for this generality is that while all the
programs given in this article are executable (and have been tested),
they run too slowly to be regarded as production simulation programs.
For example, the Poisson solver was benchmarked on some small test
cases, and compared to a straightforwardly coded Fortran equivalent.
The Haskell program ran several ($\geq 4$) orders of magnitude more slowly
than the Fortran.  This does not come as a surprise.
Similar relative performances are reported in \cite{scientific4}.  Most of the
performance degradation can be attributed to the deliberate adoption of
a very high level of abstraction in the Haskell programs, with
little or no concern for efficient execution.  (In our experience
Haskell code written in Fortran-like {\em style} is less extremely slow.
But adopting a Fortran style of programming somewhat defeats
the object of functional programming.)

     Unless compilers for functional languages make very rapid
progress, it is difficult to imagine that the kind of programs
given in this article could be competitive in the r\^{o}le of
number-crunching production codes.  The programs themselves could
be optimised to gain significant factors in speed (and efficiency of
memory utilisation).  Higher order functions could be used more
sparingly, and probably ``granularity'' could be increased in various
ways.  A few efforts were made in this direction, but these were
hindered by difficulty in understanding the behaviour (in terms of
execution speed and memory utilisation) of the lazy programs.  Programmers
who need speed are usually happy to invest some effort in optimising
their own codes, which they often do better than compilers, provided
they have a reasonably accurate idea of what goes on during execution.
It is relatively difficult to form such a mental model for a lazy
functional language; good profiling tools may make the job
easier \cite{profiling}.

     As they stand, though, the programs capture the
essence of the algorithms in a way that a production Fortran
code probably never could.  At the very least, programs like
this must have a r\^{o}le in specification or prototyping.


\appendix

%\input{appendices}

\section{\label{formLibrary}The ``form'' library}

\begin{verbatim}
module Lattice (Cell (cellEnum, cellArray,
                      cellAccum, dimension),
                Lattice (MkLattice),
                Form,
                cells, boundary,
                form, accumForm, lattice, rank,
                (!#), (//#),
                listForm, fassocs, findices, felems,
                fmap, fmapi, fzipWith) where

-- The underlying technology of lattices and their
-- forms...  This parallels the `PreludeArray'
-- section for arrays.

infixl 9 !#
infixl 9 //#

-- The `Cell' class is the analogue of the `Ix'
-- class in `PreludeCore'.

class (Eq a) => Cell a where
  cellEnum  :: Int -> [a] -> [a]
  cellArray :: Int -> [a] -> [Assoc a b] -> a -> b
  cellAccum :: Int -> [a] -> [Assoc a b] -> a -> [b]
  dimension :: a -> Int  

-- First field in `Lattice' defines extent of the
-- complex, where required. Second field is the
-- `boundary' function for the complex.

data (Cell a) =>
  Lattice a =
    MkLattice (Int -> [a]) (a -> [Assoc a Int])
data (Cell a) =>
  Form a b  = MkForm (Lattice a) Int (a -> b)

cells     :: (Cell a) => Lattice a -> Int -> [a]
boundary  :: (Cell a, Num b) =>
             Lattice a -> a -> [Assoc a b]


form      :: (Cell a) =>
   Lattice a -> Int -> [Assoc a b] -> Form a b
accumForm :: (Cell a) =>
   Lattice a -> Int -> [Assoc a b] -> Form a [b]
lattice   :: (Cell a) => Form a b -> Lattice a
rank      :: (Cell a) => Form a b -> Int
(!#)      :: (Cell a) => Form a b -> a -> b

(//#)     :: (Cell a) =>
             Form a b -> [Assoc a b] -> Form a b

listForm  :: (Cell a) =>
             (Lattice a) -> Int -> [b] -> Form a b

fassocs   :: (Cell a) => Form a b -> [Assoc a b]
findices  :: (Cell a) => Form a b -> [a]
felems    :: (Cell a) => Form a b -> [b]

fmap      :: (Cell a) => (b -> c) ->
                         Form a b -> Form a c
fmapi     :: (Cell a) => (a -> b -> c) ->
                         Form a b -> Form a c
fzipWith  :: (Cell a) => (b -> c -> d) -> Form a b ->
                         Form a c -> Form a d

-- `Lattice' primitives...

cells    (MkLattice cls _) r = cellEnum r (cls r)
boundary (MkLattice _ bdy) i =
  [j := fromIntegral s | j := s <- bdy i]
 
-- `Form' primitives...

form      l@(MkLattice cls _) r ivs =
    MkForm l r (cellArray r (cls r) ivs)
accumForm l@(MkLattice cls _) r ivs =
    MkForm l r (cellAccum r (cls r) ivs)
lattice (MkForm l _ _) = l
rank    (MkForm _ r _) = r
(!#)    (MkForm _ _ d) = d

-- Higher level operations on Forms...

a //# us   = form (lattice a) (rank a)
    ([i := a !# i
	| i <- findices a \\ [i | i := _ <- us]]
      ++ us)

listForm l r vs =
    form l r (zipWith (:=) (cells l r) vs)

fassocs  a =
    [i := a !# i | i <- cells (lattice a) (rank a)]
findices a = cells (lattice a) (rank a)
felems   a = [a !# i | i <- findices a]

fmap   f a =
    form l r [i := f (a !# i) | i <- cells l r]
    where l = lattice a
          r = rank    a

fmapi  f a =
    form l r [i := f i (a !# i) | i <- cells l r]
    where l = lattice a
          r = rank    a

fzipWith f a b =
    form l r [i := f (a !# i) (b !# i)
                | i <- cells l r]
    where l = lattice a
          r = rank    a
\end{verbatim}

\section{\label{cellInstance}An example instance of the ``Cell'' class}

\begin{verbatim}
module Cell2d (Cell2d (Cell2d0, Cell2d1, Cell2d2)) where

-- A `Cell' instance suitable for representing
-- cells in a 2d grid.

import Lattice (Cell (cellEnum, cellArray, cellAccum))

data Cell2d = Cell2d0 (Int, Int)     |
              Cell2d1 (Int, Int) Int |
              Cell2d2 (Int, Int)       deriving (Text)

ix0 (Cell2d0 i)        = i
ix1 (Cell2d1 (x, y) d) = (x, y, d)
ix2 (Cell2d2 i)        = i

bd0 [c1, c2]    = (ix0 c1, ix0 c2)
bd1 [c1, c2]    = (ix1 c1, ix1 c2)
bd2 [c1, c2]    = (ix2 c1, ix2 c2)

cell0 i         = Cell2d0 i
cell1 (x, y, d) = Cell2d1 (x,y) d
cell2 i         = Cell2d2 i

instance Cell Cell2d where

  cellEnum r cls =
    case r of
      0 -> [cell0 i | i <- range (bd0 cls)]
      1 -> [cell1 i | i <- range (bd1 cls)]
      2 -> [cell2 i | i <- range (bd2 cls)]

  cellArray r cls ivs =
    case r of
      0 -> (array (bd0 cls)
            [ix0 i := v | i := v <- ivs] !) . ix0
      1 -> (array (bd1 cls)
            [ix1 i := v | i := v <- ivs] !) . ix1
      2 -> (array (bd2 cls)
            [ix2 i := v | i := v <- ivs] !) . ix2

  cellAccum r cls ivs =
    case r of
      0 -> (accumArray (flip (:)) [] (bd0 cls)
            [ix0 i := v | i := v <- ivs] !) . ix0
      1 -> (accumArray (flip (:)) [] (bd1 cls)
            [ix1 i := v | i := v <- ivs] !) . ix1
      2 -> (accumArray (flip (:)) [] (bd2 cls)
            [ix2 i := v | i := v <- ivs] !) . ix2
\end{verbatim}

\section{\label{topologyExample}Example definition of lattice geometry}

\begin{verbatim}
module Torus2d (torus2d) where

-- `torus2d' returns an object of type `Lattice Cell2d'
-- implementing a 2 dimensional grid with
-- periodic boundary conditions.

import Lattice (Cell (cellEnum, cellArray, cellAccum),
                Lattice (MkLattice))

import Cell2d (Cell2d (Cell2d0, Cell2d1, Cell2d2))

torus2d :: Int -> Int -> Lattice Cell2d
torus2d lx ly = MkLattice (cls lx ly) (bdy lx ly)

-- Lattice bounds, specified by smallest and largest
-- cell of each dimension:

cls lx ly 0 = [Cell2d0 (0, 0),
               Cell2d0 (lx - 1, ly - 1)]
cls lx ly 1 = [Cell2d1 (0, 0) 0,
               Cell2d1 (lx - 1, ly - 1) 1]
cls lx ly 2 = [Cell2d2 (0, 0),
               Cell2d2 (lx - 1, ly - 1)]

-- Lattice topology (cell boundary):

bdy _ _ (Cell2d0 _) = []

bdy lx ly (Cell2d1 (x, y) d) =
  [Cell2d0 (x, y) :=  -1,
   Cell2d0 (case d of
              0 -> ((x + 1) `mod` lx, y)
              1 -> (x, (y + 1) `mod` ly)) := 1]

bdy lx ly (Cell2d2 (x, y)) =
  [Cell2d1 (x, y) 0 :=  1,
   Cell2d1 ((x + 1) `mod` lx, y) 1  :=  1,
   Cell2d1 (x, (y + 1) `mod` ly) 0  := -1,
   Cell2d1 (x, y) 1 := -1]
\end{verbatim}

\section{\label{running} Running some of the codes}

     In the main text detailed code was given for evolving
the states of various systems, but initialisation code was
omitted.  As an example, this appendix shows how to prepare
an initial state for the XY model simulation of section
\ref{XYimplementation}.  This is not quite trivial, because
the state incorporates an infinite list of pseudorandom
numbers at each site.  We also outline how the result of
a simulation could be used.

     Suppose $seed$ is an integer used as a random number seed, and $n$
is the size of a square lattice.  Also $theta$ is a normalised
temperature parameter.  In a physicist's simulation $theta$ might be
around unity and $n$ might be somewhere in the range 10-100.  Our
unoptimised Haskell program would certainly be limited to the bottom
end of this range.  To create a lattice with periodic boundary
conditions we define
\begin{tabbing}
\verb$  $\= \kill
\>\verb$l   = torus2d $$n$\verb$ $$n$
\end{tabbing}
Then a suitable initial state, fed to \verb$evolve$, is
\begin{tabbing}
\verb$  $\= \kill
\>\verb$ini = listForm l 0 [MkStat $$theta$\verb$ r (1.0 :+ 0.0) |$ \\
\>\verb$                    r <- ranStreams nvars $$seed$\verb$]$ \\
\>\verb$      where nvars = length (felements ini)$
\end{tabbing}
In this initial state all the $u$ variables are set to unity.
\verb$nvars$ is the number of these variables.  For a square lattice it
will just be $n^2$, but for generality it is calculated as the number
of elements of \verb$ini$.  \verb$ranStreams nvars $$seed$ returns a
list of \verb$nvars$ independent random number streams.  Our
leapfrogged linear congruential implementation is given below.  Note
the use of arbitrary precision \verb$Integer$ arithmetic.  Randoms are
scaled to the interval $[0, 1)$.
\begin{verbatim}
module Random (ranStreams) where
 
-- Random number generation.
 
a = 6364136223846793005 :: Integer
b = 1                   :: Integer
r = 2^64                :: Integer
 
next :: (Integer, Integer) -> Integer -> Integer
next (m, c) x = ((m * x) + c) `mod` r
 
intRanStream :: Int -> [Integer]
intRanStream seed = iterate (next (a, b)) (toInteger seed)
 
leap :: Int -> (Integer, Integer) -> (Integer, Integer)
-- Evaluate parameters of n-fold leapfrog generator, given
-- parameters (a, b) of simple linear congruential generator.
leap  0       _     = (1, 0)
leap (n + 1) (a, b) = ((a * m) `mod` r, (b * m + c) `mod` r)
                      where (m, c) = leap n (a, b)
 
intRanStreams :: Int -> Int -> [[Integer]]
intRanStreams n seed = map (iterate (next (leap n (a, b)))) seeds
                       where seeds = take n (intRanStream seed)
 
scale :: (Fractional a) => Integer -> a
scale n = (fromInteger n) / (fromInteger r)
 
ranStreams :: Int -> Int -> [[Float]]
ranStreams n = (map (map scale)) . (intRanStreams n)
\end{verbatim}
Typically, in this kind of simulation, the object is to find
the mean value of some ``observable'' function, averaged
over the set of random configurations generated.
For example, the {\em internal energy} is the average over
configurations of the function defined in section \ref{xy}, equation
(\ref{XYEnergy}).  This function can be computed from a form of complex
$u$ values by
\begin{verbatim}
  energy :: (Cell a) => Form a (Complex Float) -> Float
  energy u = sum [let [j := _, k := _] = boundary l i in
                    realPart ((u !# j) * conjugate (u !# k))
                 | i <- cells l 1]
             where l = lattice u
\end{verbatim}
\verb$i$ is instantiated to all links in the lattice, and the two values at
sites on the boundary of a link are multiplied together.

    Applying the map
\begin{verbatim}
  map (listForm l 0) . transpose . felems
\end{verbatim}
to the result of \verb$evolve$ produces an infinite list of arrays,
each array a complete configuration.  Applying \verb$fmap siteObs$
to each such array produces a valid argument for \verb$energy$.

     A simple-minded approach to estimating the internal energy is to
\verb$take$ a large enough finite list from the front of this stream,
preferably \verb$drop$ the first few configurations, which are
probably not sufficiently ``thermalised'', map the function that
calculates the energy of a configuration over the truncated list, and
\verb$sum$ the result (then divide by the length of the list).


\bibliography{paper}

\end{document}

