
\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.

