
\documentstyle{csart}
%\pagestyle{empty}
\sloppy
\parskip=1ex
\parindent=0em
%--------------------------------------------------------------------
% Expressing some scientific, lattice-based problems in Haskell
% -------------------------------------------------------------
% CCG seminar, 27/1/93, 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
% `dvips' went into an infinite loop...but it had produced the
% right postscript file when I killed it after a bit.
%--------------------------------------------------------------------
\setlength{\textwidth}{9cm}
\setlength{\textheight}{13cm}
\setlength{\oddsidemargin}{0cm}
\setlength{\evensidemargin}{\oddsidemargin}
\setlength{\topmargin}{-1.3cm}

\begin{document}

\title{Expressing some scientific, lattice-based problems in Haskell}
 
\author{D.B. Carpenter}

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

%\tableofcontents

\newpage

\section{Introduction}

    Will describe work done in the FAST project.
(Functional programming on Arrays of Transputers).
Aim was to produce a set of ``scientific'' demonstrator
applications.

    One of the paradigms easily expressible in lazy functional languages
is process-parallel programming---FAST project produced
a system for parallel execution of annotated programs
written in this style.

\subsection*{Plan of talk}
\begin{enumerate}
\item
Introduction to functional programming and Haskell
\item
Lattices and generalised arrays
\item
Implementation of some scientific algorithms
\end{enumerate}

\newpage

\section{About functional languages}

     Functional programming languages are very mathematically
``pure''.  They have no intrinsic idea of control-flow, store,
assignment, etc.  A functional program is a mathematical expression.
This expression may be a function---say from inputs to outputs.

     Higher order functions (functions which take other functions as
parameters, or return other functions as results) can be used freely.
In practise this allows imperative control-flow to be emulated, if
necessary.

     Lazy evaluation permits manipulation of ``infinite'' data
structures, in particular infinite lists.  These provide an alternate
means of representing state transition---the history of a transition
system is an infinite list (stream) of system states.

     So we can recover sequential programming if desired (usually
not), and also model more general imperative paradigms, including
process parallelism, easily.

     Beyond that, functional language are very {\em uniform}, and free
of ad hoc restrictions.  This makes them very expressive, and allows
high levels of abstraction.

\newpage

\section{Some basic types}

     The Haskell system comes with a large number of predefined types
(including \verb$Int$, \verb$Integer$, \verb$Bool$, \verb$Char$,
\verb$Float$, complex types, etc).  A constant can be specified
to have a particular type by a declaration
\begin{verbatim}
  a :: Int
\end{verbatim}

     A function type is a type like any other
\begin{verbatim}
  f :: Int -> Float
\end{verbatim}
means \verb$f$ a function from \verb$Int$ to \verb$Float$.
Functions can take functions as arguments, or return functions
as results
\begin{verbatim}
  g :: (Int -> Float) -> Float
  h :: Int -> (Int -> Float)
  h :: Int -> Int -> Float
\end{verbatim}
Multiple arguments represented by ``currying''.
The two signatures of \verb$h$ are equivalent by definition
and represent a function from an \verb$Int$ to a \verb$Int -> Float$---aka
a function from two \verb$Int$s to a \verb$Float$.

     Another basic type is the tuple
\begin{verbatim}
  a :: (Int, Float)
  b :: (Int, Int, Int)
\end{verbatim}
A function from a pair to a float:
\begin{verbatim}
  h' :: (Int, Int) -> Float
\end{verbatim}

\newpage

\section{Function application---expressions}

    Function application represented by juxtaposition
\begin{verbatim}
  f 3
  g f
  h 3 4
  h' (3, 4)
\end{verbatim}
Functions of two arguments can be used as ``infix'' operators
\begin{verbatim}
  (+), (*) :: Float -> Float -> Float
\end{verbatim}
(obviously these are predefined and overloaded in practise)
so arithmetic expressions can be written in familiar ways
(and other function applications can be written in unfamiliar ones)
\begin{verbatim}
  3.0 + (4.0 * f 3)
  3 `h` 4
\end{verbatim}
Due to currying, functions can be partially applied
\begin{verbatim}
  h 3
\end{verbatim}
is a function from \verb$Int$ to \verb$Float$, and
\begin{verbatim}
  (2.0 *)
\end{verbatim}
is a function from \verb$Float$ to \verb$Float$.  (These
are elementary examples of function-valued expressions).

\section{Data constructor function definition}

     ``Algebraic types'' are roughly equivalent to
``variant records'' or ``Unions'' in some languages.
They have an associated set of ``constructor functions'',
which are genuine functions, but with limited power.
\begin{verbatim}
  data IntList = Node Int IntList | Nil
\end{verbatim}
defines an algebraic type \verb$IntList$ with
two constructor functions \verb$Node$ and \verb$Nil$.
\verb$Node$ is of type \verb$Int -> IntList -> IntList$.
\verb$Nil$ is a constant (``function of zero arguments'')
of type \verb$IntList$.

     The following expressions are of type \verb$IntList$
\begin{verbatim}
  Nil
  Node 1 (Node 2 (Node 3 Nil)))
  1 `Node` Nil
\end{verbatim}

     Not surprisingly, Haskell has built in syntax for representing
general lists---this is just given as an example of algebraic types.

\newpage

\section{General function definition}

     General functions are defined by a straightforward syntax
\begin{verbatim}
  double :: Float -> Float
  double x = 2.0 * x

  factorial :: Int -> Int
  factorial n = if n == 0
                  then 1
                  else n * factorial (n - 1)
\end{verbatim}
Functions on algebraic types are handled by ``pattern-matching''
\begin{verbatim}
  isEmpty :: IntList -> Bool
  isEmpty (Node x xs) = False
  isEmpty  Nil        = True

  car :: IntList -> Int       -- ``head''
  car (Node x xs) = x

  cdr :: IntList -> Intlist   -- ``tail''
  cdr (Node x xs) = xs
\end{verbatim}
Note \verb$car$ and \verb$cdr$ are undefined if the argument is
\verb$Nil$.

     \verb$double$ could also be defined by
\begin{verbatim}
  double = (2.0 *)
  double = (\x -> 2.0 * x)
\end{verbatim}

\newpage

\section{Special syntax for lists}

     Lists can be defined from scratch---a restricted
form (\verb$Int$s only) was already given.  They are
used so commonly that some special syntax has been added
to the language.  The equivalent of \verb$Node$ is \verb$(:)$,
and it's normally used as an infix function.  The equivalent
of \verb$Nil$ is \verb$[]$.  Additionally a list of several
items can be written as comma-separated sequences in square
brackets.  The following are equivalent:
\begin{verbatim}
  (:) 1.0 ((:) 2.0 [])
  1.0 : (2.0 : [])
  [1.0, 2.0]
\end{verbatim}
This is a list of two \verb$Float$ items.  Its {\em type} is
written: \verb$[Float]$.  This is an example of a parametrised
type---in this case the parameter type is \verb$Float$.

     A ``list comprehension'' is another special piece of
syntax for building lists.  A simple example is the expression
\begin{verbatim}
  [x^2 | x <- xs]
\end{verbatim}
Here \verb$xs$ is a list of numeric values.  The result of
the expression is a list formed by instantiating \verb$x$
to the value of each item in \verb$xs$, and computing the
expression preceding the \verb$|$.
If \verb$xs = [1, 2, 3]$, the result of the above expression is
\verb$[1, 4, 9]$.  (Fortran users might interpret this as something
like an implied do loop!)  This syntax has a few extensions---multiple
generators can be placed after the \verb$|$ (nested implied
do loop!), and predicates can be added to filter out certain
values of the ``dummy indices''.

\newpage

\section{Normal order evaluation}

     Lazy languages implement ``normal order evaluation'' of
expressions.  Evaluation of an expression starts ``at the
outermost level'', and no sub-expression is evaluated
unless it is actually needed. Eg, if we had
\begin{verbatim}
  cond :: Bool -> a -> a -> a
  cond True  x y = x
  cond False x y = y

  n = 0
\end{verbatim}
then the expression \verb$cond (n /= 0) (1 / n) 0$ is well-defined.
Only the first argument has to be evaluated to make the pattern-match.

\newpage

\section{Infinite data structures}

     This means that we can deal in expressions which would
take $\infty$ly long to evaluate in full, so long as we only
choose to look at a parts of the result which can be evaluated
in finite time.  It often simplifies construction of data structures
if one doesn't worry initially about termination.

     An application of infinite lists is in modelling
state transition.  Define
\begin{verbatim}
  process :: (s -> e -> s) -> s -> [e] -> [s]
  process trans st ~(env:envs) =
        st : process trans (trans st env) envs
\end{verbatim}
This takes: a transition function from a state and a current
``environment value'' to a new state; an initial state; and an $\infty$
list of successive environment values.  It returns an $\infty$
list of states.

    As an example, take the ``transition function'' to be \verb$(+)$.
This adds the current value of the ``environment'' to the current value
of the ``state''.  Then we would have
\begin{verbatim}
  process (+) 0 [1, 3, 5, 7, ...]  ==
                        [0, 1, 4, 9, 16, ...]
\end{verbatim}
So long as the computation only ultimately requires a finite number of
items from this list, the fact that it formally
contains an $\infty$ number of items doesn't matter.

\section{Some standard functions}

     Some standard functions:
\begin{verbatim}
  map f xs
     = [f x | x <- xs]

  zipWith f [x, y, ...] [a, b, ...]
     = [f x a, f y b, ...]

  transpose [[a1, a2, ...], [b1, b2, ...], ...]
     = [[a1, b1, ...], [a2, b2, ...], ...]

  (f . g) x
     = f (g x)
\end{verbatim}
All these can easily be defined in terms of the operations
given so far (and are so-defined in the prelude).
Where they take lists, all of them are applicable to finite or infinite lists.

\newpage

\section{Haskell arrays}

     An array can be constructed from an ``association list''.
An association is a pair constructed with the
infix data constructor \verb$(:=)$, binding a subscript to
a value.  \verb$array$ is a
function which takes a pair of bounds and an association
list, and returns an array.  Eg
\begin{verbatim}
  a = array (1, 3) [1 := Red,
                    2 := Blue, 3 := Green]
\end{verbatim}
is an array with subscript type \verb$Int$ (tuples are
also allowed, providing multi-dimensional arrays)
and values in some algebraic type representing colours.
The subscripting operator is the infix function \verb$(!)$, so
\begin{verbatim}
  a ! 2  ==  Blue
\end{verbatim}
Arrays aren't ``primitive'' in Haskell.  They can be implemented
(inefficiently) in terms of function.  They don't even have any special
syntax, as lists do.  So it's possible to implement ones own array
system (provided efficiency isn't a overiding concern).

\newpage

\section{Lattices}

     QCD nomenclature for the ``cells'' of a lattice
\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}
}
Variables in gauge theory and other lattice-based
models ``live'' on different dimensional cells of
the lattice.

     Similar ``cell complexes'' arise in topology.  The
structure of a complex can be encoded in a ``boundary
operator'', $\partial$.
\centerline{
%\begin{picture}(220, 130)(-20, -30)
\begin{picture}(120, 110)(-20, -10)
\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}
}
\begin{eqnarray*}
\partial F & = & G - H \\
\partial A & = & B + C - D - E \\
\partial H & = & 0
\end{eqnarray*}
The boundary operation is powerful enough to encode the topology
of an {\em arbitrary} cell complex.

\newpage

\section{Cells and generalised arrays in Haskell}

Cells of a 2d square lattice can be represented by
\begin{verbatim}
  data Cell2d = Cell2d0 (Int, Int)     |
                Cell2d1 (Int, Int) Int |
                Cell2d2 (Int, Int)
\end{verbatim}
A ``form'' is an array whose index space is the set of
cells of a particular dimension in a particular lattice.

Suppose \verb$l$ is a data structure representing a $2 \times 2$
lattice
\begin{verbatim}
  a = form l 0 [Cell2d0 (0, 0) := Red,
                Cell2d0 (0, 1) := Blue,
                Cell2d0 (1, 0) := Green,
                Cell2d0 (1, 1) := Yellow]
\end{verbatim}
is a ``$0$-form'', taking values (colours, here) on the
{\em sites} of \verb$l$.  Subscript by \verb$(!#)$:
\begin{verbatim}
  a !# Cell2d0 (1, 0)  ==  Green
\end{verbatim}

Haskell has a class system, which allows the basic ``form''
operations to be {\em overloaded} over arbitrary instances
of ``cell'' types and so arbitrary lattice shapes---\verb$Cell2d$
is just one example.  The type signature of \verb$form$ is
\begin{verbatim}
  form :: (Cell a) =>
    Lattice a -> Int -> [Assoc a b] -> Form a b
\end{verbatim}

\section{More operations on cells and forms}

\begin{itemize}
\item{\verb$cells l r$}
returns a list of all cells of dimension
\verb$r$ in lattice \verb$l$.

\item{\verb$boundary l i$}
returns a list of the cells (paired
with their incidence function) on the boundary of cell \verb$i$
in lattice \verb$l$.

\item{\verb$fmap$}
is like \verb$map$.  It takes a mapping function
and a form, and returns a form.

\item{\verb$fzipWith$}
is like \verb$zipWith$.  It takes a zipping function
and a pair of forms, and returns a form.

\item{\verb$nns$}
is complicated.
\end{itemize}

\verb$nns$ is a function which takes a mapping function and
a $0$-form, and returns a $0$-form in which each element is
a mapped list of the nearest neighbour elements of the
original form.  If $j_1, j_2, ...$ are the neighbours of
site $i$, and \verb$b = nns f a$, then
\newline
\newline
\verb$  b !# $$i$\verb$ == [f $$i$\verb$ $$j_1$\verb$ (a !# $$j_1$\verb$), f $$i$\verb$ $$j_2$\verb$ (a !# $$j_2$\verb$), ...]$
\newline


\newpage

\section{Implementation of {\bf nns}}

     All operations are ``local'' except \verb$boundary$.
This works for the sites of {\em any} shape lattice.
\begin{verbatim}
  mltply m ivs = [i := m * v | i := v <- ivs]

  neighbours l r = accumForm l r
    [j := mltply v
       (if v == (-1) then reverse b else b)
      | b      <- [boundary l i |
                   i <- cells l (r + 1)],
        j := v <- b]

  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}

\newpage

\section{Poisson relaxation}

Continuum Poisson equation:
\begin{displaymath}
(\nabla^2) (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.

\newpage

A type to capture the state of a site
\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}
Given an array containing streams at each site,
representing the history of that site, \verb$env$
returns an array of streams of each site's environment.
\begin{verbatim}
env = fmap transpose . nns (\_ _ -> map siteObs)
\end{verbatim}
Given an array of initial states, evolve produces
an array of site ``history'' streams
\begin{verbatim}
evolve :: Form Cell2d State -> Form Cell2d [State]
evolve init = x
 where x = fzipWith (process siteTrans) init e
       e = env x
\end{verbatim}
 
\newpage

\section{Red-black relaxation}

\begin{verbatim}
colour (Cell2d0 (x, y)) = (x + y) `mod` 2
 
env = fmap transpose .
        nns (\i j -> map siteObs .
                     if colour j < colour i
                       then tail
                       else id)
\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 ``vortices''.

\newpage

Only changes to Poisson code are in definition of {\em State} type.

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

\section{Lattice gauge theory}

     The difference between $U(1)$ gauge theory (lattice QED) and the XY
model is geometric.  In gauge theory variables associated with links
of the lattice, and the energy function is computed on the
plaquettes. 
\centerline{
\begin{picture}(120, 110)(0, -25)
\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}
}
Associate a group element
$U_P$ with any path $P$ in the lattice by multiplying
together the link variables along the path.  Energy function is then
\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.

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

\newpage

     Replace \verb$env$ by
\begin{verbatim}
  staples f a =
    fmapi (\i ivss ->
           [f i [(j, v, a !# j) | j := v <- snip i ivs]
           | ivs <- ivss]) (neighbours l r)
    where l           = lattice a
          r           = rank a
          snip i ivs  =
            tail (rotate (\(j := _) -> j == i) ivs)
          rotate p xs = b ++ f where (f, b) = break p xs

  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}
and \verb$colour$ and \verb$deltaE$ by
\begin{verbatim}
  colour (Cell4d1 (t, x, y, z) d) =
      4 * ((t + x + y + z) `mod` 2) + d.

  deltaE neighbs s s' =
    - real (sum neighbs * (s' - s))
      where real (x :+ _) = x
\end{verbatim}


\end{document}
