
\documentstyle{csart}
%\pagestyle{empty}
\sloppy
\parskip=1ex
\parindent=0em
%--------------------------------------------------------------------
% Expressing some grid-based algorithms in Haskell
% ------------------------------------------------
% FAST meeting, 30/11/92, Southampton.
% Latex format borrowed from JHM.
%
% To obtain slides (portrait format) do:
%   latex progs
%   dvips -x 1800 -o slides.ps slides.dvi
%   lpr slides.ps
%--------------------------------------------------------------------
\setlength{\textwidth}{9cm}
\setlength{\textheight}{13cm}
\setlength{\oddsidemargin}{0cm}
\setlength{\evensidemargin}{\oddsidemargin}
\setlength{\topmargin}{-1.3cm}

\begin{document}

\title{Expressing some grid-based algorithms in Haskell}
 
\author{D.B. Carpenter}

%\address{Department of Electronics and Computer Science, \\
%        Southampton University, \\
%         Southampton, S09 5NH} 
        
\maketitle

\tableofcontents

\newpage

\section{Ojectives}

\begin{itemize}
\item
Explore expressive power of functional language,
applied to some favourite scientific problems.
\item
Exercise the FAST compiler, on some number-crunching
applications.
\item
Produce some demos for Hugh.
\end{itemize}

The applications will be grid oriented, including
PDE (Poisson equation), statistical mechanics spin
models, and lattice gauge theory (QCD).

First aim is to exploit the high-level features of
functional language to express the algorithms in
the most natural mathematical way.

Second aim is to go fast.

\newpage

\section{Nomenclature}

Will use nomenclature common in QCD:

\begin{description}

\item{\bf Lattice} = ``grid''

\begin{picture}(40, 40)(-10, -10)
\thicklines
\multiput(-5,  0)( 0, 10){3}{\line(1, 0){30}}
\multiput( 0, -5)(10,  0){3}{\line(0, 1){30}}
\end{picture}

\item{\bf Site} = ``point'' or ``node''.

\begin{picture}(40, 40)(-10, -10)
\multiput(-5,  0)( 0, 10){3}{\line(1, 0){30}}
\multiput( 0, -5)(10,  0){3}{\line(0, 1){30}}
\put(0, 0){\circle*{4}}
\end{picture}

\item{\bf Link} = ``edge'' or ``arc''.

\begin{picture}(40, 40)(-10, -10)
\multiput(-5,  0)( 0, 10){3}{\line(1, 0){30}}
\multiput( 0, -5)(10,  0){3}{\line(0, 1){30}}
\thicklines
\put(0,  0){\line(1, 0){10}}
\end{picture}

\item{\bf Plaquette} = ``square''.

\begin{picture}(40, 40)(-10, -10)
\multiput(-5,  0)( 0, 10){3}{\line(1, 0){30}}
\multiput( 0, -5)(10,  0){3}{\line(0, 1){30}}
\thicklines
\multiput( 0,  0)( 0, 10){2}{\line(1, 0){10}}
\multiput( 0,  0)(10,  0){2}{\line(0, 1){10}}
\end{picture}

\end{description}

In general, lattice has arbitrary dimensionality, arbitrary
boundary-conditions (eg, open or toroidal).

Can regard lattice simply as an array of sites,
or as some sort of complex of sites, links, plaquettes,
etc, with some incidence function between these elements.

\section{Poisson relaxation}

Continuum Poisson equation:
\begin{displaymath}
(\Box u) (x, y) = c (x, y),
\end{displaymath}
where $c$ is given charge distribution; $u$ is electrostatic potential,
to be calculated.  Discretise:
\begin{displaymath}
\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{displaymath}
Naive relaxation---rearrange equation:
\begin{displaymath}
u (x, y) = \frac{1}{4} (\mbox{neighbouring $u$ values} - c (x, y)),
\end{displaymath}
then iterate from some initial guess---new $u$ obtained by substituting
old $u$ in RHS.

    Two natural ways of arranging this iteration in a lazy language:
\begin{itemize}
\item
Build a stream of arrays.  Each successive item in the stream is
an improved approximation to the array $u$.
\item
Build an array of streams.  Each element of the array is a list
of gradually improving approximations to the value of $u$ at
a particular site.
\end{itemize}

\newpage

Introduce some useful bits and pieces.  First
some very general operations:

A handy generalisation of {\em amap}
\begin{verbatim}
  amapi f a = array b [i := f i (a ! i) |
                       i <- range b]
              where b = bounds a
\end{verbatim}

A useful primitive for evolving state transition systems
\begin{verbatim}
  process :: (s -> e -> s) -> s -> [e] -> [s]
  process trans st (env:envs) =
        st : process trans (trans st env) envs
\end{verbatim}
(NB: another implementation is {\em process} = {\em scanl})

\newpage

Some things more specific to Poisson's equation:

An operation for accessing neighbours
\begin{verbatim}
  neighbours (x, y) = [(x - 1, y), (x + 1, y),
                       (x, y - 1), (x, y + 1)]

  a !# i = [a ! j | j <- neighbours i]
\end{verbatim}
A type to capture the state of a site
\begin{verbatim}
  data State = Interior Float Float |
               Boundary       Float

  siteObs (Interior _ u) = u
  siteObs (Boundary   u) = u

  siteTrans (Interior c _) nbs =
             Interior c (0.25 * (sum nbs - c))
  siteTrans b              _   = b
\end{verbatim}
{\em siteObs} extracts the potential.
{\em siteTrans} is the transition function for
relaxation (sites on a boundary aren't updated).

\newpage

A function which returns the history of a site,
given an initial state and the history of its neighbours:
\begin{verbatim}
  siteStrm _ ini nbStrms =
                process siteTrans ini env
                where env = transpose nbStrms
\end{verbatim}
Now tie up the site streams:
\begin{verbatim}
  arrStrms inia = ssa
      where ssa     = amapi f inia
            f i ini = siteStrm i ini (osa !# i)
            osa     = amap (map siteObs) ssa
\end{verbatim}
(Use convention that identifiers ending $s$ are lists, in $a$ are
arrays; $s$ stands for {\em State}, $o$ for {\em observable}).
{\em arrStrms} returns an array of {\em State}
streams, given an array of initial states.

Note---neighbour operations don't go out of array bounds,
due to lazy evaluation.

For comparision, a ``stream of arrays'' implementation built from
the same components:
\begin{verbatim}
  arrTrans sa oa = amapi f sa
           where f i s = siteTrans s (oa !# i)

  strmArrs inia = sas
        where sas   = process arrTrans inia oas
              oas   = map (amap siteObs) sas
\end{verbatim}

\newpage

\section{XY model}

     Solving Poisson equation is equivalent to the variational
problem: minimise
\begin{displaymath}
  E (u) = \sum_{x, y} \frac{1}{2} (
  \begin{array}[t]{l}
  (u (x + 1, y) - u (x, y))^2 + \\
  (u (x, y + 1) - u (x, y))^2 \,\,) + \\
  c (x, y) u (x, y).
  \end{array}
\end{displaymath}
Typical statistical physics simulations involves a similar
type of energy
functional, but rather than minimising $E$, the object is
to generate an ensemble of configurations with a Boltzmann
probability distribution:
\begin{displaymath}
  P (u) \propto \exp (- E (u) / \theta)
\end{displaymath}
where $\theta$ is the temperature of the system.
(NB.  Equivalent to minimisation when
$\theta \rightarrow 0$; Simulated Annealing).

     XY model is a ``spin system''.  The degree of freedom
at a site is a unit 2-vector, or unit magnitude complex number
(rather than a single real number in Poisson).  Energy function is
\begin{displaymath}
  E (u) = - \Re \sum_{x, y} u^* (x, y) u (x + 1, y) +
                            u^* (x, y) u (x, y + 1)
\end{displaymath}

     Model is famous for having an unusual phase transition,
driven by condensation of topological dislocations
(vortices).

\newpage

Main changes to Poisson code are in definition of {\em State} type.
(I also changed \verb$!#$ to implement periodic boundary
conditions, so no worry about boundaries and interiors here.)

Initial array has to include a stream of independent random
numbers for each {\em State} variable.

\begin{verbatim}
  data State = MkStat Float [Float] Spin 
  -- temperature, a list of randoms and ``spin''.

  type Spin  = Complex Float

  siteObs :: State -> Spin
  siteObs (MkStat _ _ s) = s

  siteTrans :: State -> [Spin] -> State
  -- Local Metropolis update
  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}

\newpage

Details of Metropolis code:
\begin{verbatim}
  -- Energy difference between two spin states.
  deltaE neighbs s s' = x
    where x :+ _ = - (conjugate (sum neighbs)) *
                     (s' - s)

  -- Canditate 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}

\newpage

Would like to automatically spot vortex configurations, eg

\begin{picture}(120, 100)(-20, -20)
\multiput(-10,   0)( 0, 20){4}{\line(1, 0){80}}
\multiput(  0, -10)(20,  0){4}{\line(0, 1){80}}
\thicklines
\put( 0,  0){\line(-1, 1){10}}
\put( 0, 20){\line(-1, 2){10}}
\put( 0, 40){\line( 1, 2){10}}
\put( 0, 60){\line( 1, 1){10}}
\put(20,  0){\line(-2, 1){10}}
\put(20, 20){\line(-1, 1){10}}
\put(20, 40){\line( 1, 1){10}}
\put(20, 60){\line( 2, 1){10}}
\put(40,  0){\line(-2,-1){10}}
\put(40, 20){\line(-1,-1){10}}
\put(40, 40){\line( 1,-1){10}}
\put(40, 60){\line( 2,-1){10}}
\put(60,  0){\line(-1,-1){10}}
\put(60, 20){\line(-1,-2){10}}
\put(60, 40){\line( 1,-2){10}}
\put(60, 60){\line( 1,-1){10}}
\end{picture}
\begin{picture}(120, 100)(-20, -20)
\multiput(-10,   0)( 0, 20){4}{\line(1, 0){80}}
\multiput(  0, -10)(20,  0){4}{\line(0, 1){80}}
\thicklines
\put( 0,  0){\line(-1,-1){10}}
\put( 0, 20){\line(-2,-1){10}}
\put( 0, 40){\line(-2, 1){10}}
\put( 0, 60){\line(-1, 1){10}}
\put(20,  0){\line(-1,-2){10}}
\put(20, 20){\line(-1,-1){10}}
\put(20, 40){\line(-1, 1){10}}
\put(20, 60){\line(-1, 2){10}}
\put(40,  0){\line( 1,-2){10}}
\put(40, 20){\line( 1,-1){10}}
\put(40, 40){\line( 1, 1){10}}
\put(40, 60){\line( 1, 2){10}}
\put(60,  0){\line( 1,-1){10}}
\put(60, 20){\line( 2,-1){10}}
\put(60, 40){\line( 2, 1){10}}
\put(60, 60){\line( 1, 1){10}}
\end{picture}\\
Number of vorticies inside a given cyclic
contour given by integrating the local
rate of change of vector direction around
the contour (and dividing by $2 \pi$).

The smallest cyclic ``contour'' in lattice
is plaquette.  Number of vorticies in
plaquette:

\begin{picture}(100,  80)(-40, -20)
\multiput(-10,   0)( 0, 40){2}{\line(1, 0){60}}
\multiput(  0, -10)(40,  0){2}{\line(0, 1){60}}
\thicklines
\put( 0,  0){\line(-1, 1){20}}
\put(-10,20){$\theta_1$}
\put( 0, 40){\line( 1, 1){20}}
\put(15,45){$\theta_2$}
\put(40,  0){\line(-1,-1){20}}
\put(45,15){$\theta_3$}
\put(40, 40){\line( 1,-1){20}}
\put(15,-10){$\theta_4$}
\end{picture}\\
is
\begin{displaymath}
(\Delta (\theta_2, \theta_1) +
\Delta (\theta_3, \theta_2) +
\Delta (\theta_4, \theta_3) +
\Delta (\theta_1, \theta_4)) / 2 \pi
\end{displaymath}
where
\begin{displaymath}
\Delta(x, y) = x - y \bmod 2 \pi
\end{displaymath}
Can't be done conveniently with the simple
{\em amapi}, \verb$!#$ primitives introduced so far.

\newpage

\section{Lattice gauge theory}

Quantum field theories like QCD can be reformulated to
look like statistical mechanics problems.

In a gauge theory, the variables are {\em group elements}, $U$,
which lie on {\em links} of lattice.  Energy function
associated with plaquettes of lattice
(cf XY model, where variables associated with
sites and energy associated with links).

With a closed path $P$ starting at any point
in lattice can associate a group element $U_P$ by multiplying
together the link variables in path order.

Energy is
\begin{displaymath}
E (U) = \sum_{P \in \mbox{\em plaquettes}} \chi (U_P)
\end{displaymath}
where {\em plaquettes} is the set of plaquette-paths in the lattice
and $\chi$ is a characteristic function for the
group.  Normally $U$ represented by complex
matrices (eg, $3 \times 3$ unitary matrices) and $\chi$
just real part of trace of a matrix.

Gauge theories ``trivial'' in 2 dimensions, theoretically interesting
in 3 dimensions, and real QCD lives in 4 dimensions.

Again, {\em amapi}, \verb$!#$ primitives inadequate to
capture this geometry.

\newpage

\section{A more general view of lattice geometry?}

     Currently, effectively just have arrays defined
on sites of lattice.

\begin{itemize}
\item{Program:}
     
     Generalise so that array index
space can also be links, plaquettes, etc.

     Replace {\em neighbours} operation by operations
which take the {\em boundary} or {\em coboundary} of
a particular cell


\begin{picture}(60, 20)(-10, -10)
\put(5, -5){cell}
\end{picture}
\begin{picture}(60, 20)(-10, -10)
\put(5, -5){bnd}
\end{picture}
\begin{picture}(60, 20)(-10, -10)
\put(5, -5){cob}
\end{picture}

\begin{picture}(60, 40)(-10, -10)
\multiput(-5,  0)( 0, 10){3}{\line(1, 0){30}}
\multiput( 0, -5)(10,  0){3}{\line(0, 1){30}}
\put(10, 10){\circle*{4}}
\end{picture}
\begin{picture}(60, 40)(-10, -10)
\multiput(-5,  0)( 0, 10){3}{\line(1, 0){30}}
\multiput( 0, -5)(10,  0){3}{\line(0, 1){30}}
\end{picture}
\begin{picture}(60, 40)(-10, -10)
\multiput(-5,  0)( 0, 10){3}{\line(1, 0){30}}
\multiput( 0, -5)(10,  0){3}{\line(0, 1){30}}
\thicklines
\put(10,  10){\line(1, 0){10}}
\put(10,  11){\line(1, 0){10}}
\put(10,  10){\line(0, 1){10}}
\put( 9,  10){\line(0, 1){10}}
\put(10,  10){\line(-1, 0){10}}
\put(10,  11){\line(-1, 0){10}}
\put(10,  10){\line(0, -1){10}}
\put( 9,  10){\line(0, -1){10}}
\end{picture}

\begin{picture}(60, 40)(-10, -10)
\multiput(-5,  0)( 0, 10){3}{\line(1, 0){30}}
\multiput( 0, -5)(10,  0){3}{\line(0, 1){30}}
\thicklines
\put(10,  10){\line(0, 1){10}}
\put( 9,  10){\line(0, 1){10}}
\end{picture}
\begin{picture}(60, 40)(-10, -10)
\multiput(-5,  0)( 0, 10){3}{\line(1, 0){30}}
\multiput( 0, -5)(10,  0){3}{\line(0, 1){30}}
\put(10, 10){\circle*{4}}
\put(10, 20){\circle*{4}}
\end{picture}
\begin{picture}(60, 40)(-10, -10)
\multiput(-5,  0)( 0, 10){3}{\line(1, 0){30}}
\multiput( 0, -5)(10,  0){3}{\line(0, 1){30}}
\thicklines
\put( 2, 11){\framebox(6, 7){}}
\put(12, 11){\framebox(6, 7){}}
\end{picture}

\begin{picture}(60, 40)(-10, -10)
\multiput(-5,  0)( 0, 10){3}{\line(1, 0){30}}
\multiput( 0, -5)(10,  0){3}{\line(0, 1){30}}
\thicklines
\put(12, 11){\framebox(6, 7){}}
\end{picture}
\begin{picture}(60, 40)(-10, -10)
\multiput(-5,  0)( 0, 10){3}{\line(1, 0){30}}
\multiput( 0, -5)(10,  0){3}{\line(0, 1){30}}
\thicklines
\put(10,  10){\line(1, 0){10}}
\put(10,  11){\line(1, 0){10}}
\put(10,  10){\line(0, 1){10}}
\put( 9,  10){\line(0, 1){10}}
\put(20,  20){\line(-1, 0){10}}
\put(20,  21){\line(-1, 0){10}}
\put(20,  20){\line(0, -1){10}}
\put(19,  20){\line(0, -1){10}}
\end{picture}
\begin{picture}(60, 40)(-10, -10)
\multiput(-5,  0)( 0, 10){3}{\line(1, 0){30}}
\multiput( 0, -5)(10,  0){3}{\line(0, 1){30}}
\end{picture}

Now replace {\em amapi} by a function {\em bndmapi} which takes
an array over (eg) the plaquettes to an array over the links,
applying the map function on a ``per-incidence'' (overlap of
plaquette and link) basis.  Similarly {\em cobmapi}.

\end{itemize}

\end{document}

