\documentclass{article}
\usepackage{sober,url,amsmath}
\newcommand{\note}[1]{\emph{[#1]}}
\newcommand{\emonk}{{\sf electric monk}}
\newcommand{\ema}{{\sf em}}
\newcommand{\fixme}[1]{{\bf #1}}
\newcommand{\code}[1]{{\tt #1}}
\title{\emonk: a system for auto-deriving
spatial moment equations}
\author{Ben Bolker}
\date{\today}
\begin{document}
\newcommand{\mma}{Mathematica}
\newcommand{\xtc}{{\tt xtc}}

\maketitle
\section{Introduction}
\emonk\ (\ema\ for short: named after the character in a Douglas Adams
novel)\footnote{
``The Electric Monk was a labour-saving device, like a
dishwasher or a video recorder.  Dishwashers washed tedious dishes for
you, thus saving you the bother of washing them yourself, video
recorders watched tedious television for you, thus saving you the
bother of looking at it yourself; Electric Monks believed things for
you, thus saving you what was becoming an increasingly onerous task,
that of believing all the things the world expected you to believe'':
Douglas Adams, \emph{Dirk Gently's Holistic Detective Agency}}
is a collection of \mma\ code for automatically deriving,
analyzing, and numerically integrating spatial moment equations.
\ema's core consists of five sets of rules:
\begin{itemize}
\item{rules for finding the expected changes of variables and products of
variables, given a set of transitions and transition rates;}
\item{rules for simplifying integrals involving \emph{kernels}
(probability distributions of interaction strengths as a function of
distance), means, covariances, etc.;}
\item{rules for making \emph{moment closure} approximations
of these equations;}
\item{rules for converting the resulting moment equations from
\mma\ to input files for \xtc, a numerical integration program by Bard
Ermentrout;}
\item{rules for simplifying combinations of Fourier-transformed
Laplacian/Bessel kernels.}
\end{itemize}
This package is intended to make the process of exploring spatial
moment equations less painful.  As with any free package, it comes
with NO WARRANTY, expressed or implied.  (It is released under the Gnu
Public License, in case that means anything to you.)  As with any
software package for doing ``analysis'', it is your responsibility not
to do anything silly, and to understand what your results mean \ldots
\section{Starting up}
Start \mma\ (either in text or notebook mode) and
load in the package:
\begin{equation}
\verb|<<emonk.ma|
\end{equation}
To get help on a particular function \verb|foo| (provided I've
written it) type \verb+?foo+; for example,
\begin{equation}
\verb|?meaneq|
\end{equation}

\section{Defining a model}
Enter the transition rates for an individual-based model as a nested
list in \mma\ form as follows:
\begin{equation}
\verb+model = {rule1,rule2,rule3}+
\end{equation}
where each of the rules specifies a transition rule as defined
below.  (Note that each rule is itself a list in \mma, so
it will have its own set of curly brackets \{\} around it.)

Each transition rule is
defined as a list of
three components.  The first component is the list of
variables (states) and locations changed by the transition; the second
is the (list of) sizes of the changes to each variable made by the transition
(typically $\pm$ 1); and the third is the rate of the transition
(a formula, \emph{not} a list)
expressed as a function of local and nearby state values.

\subsection{Local (density-independent) removal/addition}
For example, a simple mortality term 
\begin{equation}
N(x) \to N(x)-1 \mbox{ at rate } \mu N(x)
\end{equation}
would be expressed in the \mma\ model as
\begin{equation}
\verb+rule1={{N[x]},{-1},{mu*N[x]}}+.
\end{equation}

\subsection{Local (density-independent) transformation}
A density-independent transformation such as germination
(from fruit (\verb+N+) to seedling (\verb+S+))
would be expressed mathematically
as 
\begin{equation}
\{N(x),S(x)\} \to \{N(x)-1,S(x)+1\} \mbox{ at rate } g N(x)
\end{equation}
and in \mma\ as
\begin{equation}
\verb|rule2={{N[x],S[x]},{-1,+1},g*N[x]}|.
\end{equation}

\subsection{Neighbourhood- or (density-)dependent addition or removal}
Transitions can also depend on neighborhood densities,
expressed as integrals over kernels.
Fruit deposition from neighboring trees (\verb+T+)
which have (passive) dispersal kernel $D$ and seed production rate $s$
would lead to a transition
\begin{equation}
N(x) \to N(x)+1 \mbox{ at rate } \int_\Omega s D(|y-x|) T(y) \, dy
\end{equation}
(where $\Omega$ is the entire space, $x$ and $y$ are vectors,
and $|x-y|$ is the distance between them)
or
\begin{equation}
\verb|rule3={{N[x]},{+1},int[s*Dk[y-x]*T[y],y]}|
\end{equation}
(the integration variable, $y$ in this case, is
given as the second argument to the integral function \verb|int|).
(Note that the kernel is also expressed as a function (using square
brackets in \mma); we use \verb|Dk| instead of \verb|D| because
the latter is a reserved word in \mma.)

\subsection{Neighbourhood- (or density-)dependent transformation}
Infection of a susceptible individual (\verb+N+)
by propagules from
infected individuals (\verb+I+) in the neighborhood is a transformation
that depends on local neighborhood density of infectives

\begin{equation}
N(x) \to N(x)-1,
I(x) \to I(x)+1,
 \mbox{ at rate } \beta N(x) \int_\Omega U(|y-x|) I(y) \, dy
\end{equation}
or
\begin{equation}
\verb|rule4={{N[x],I[x]},{-1,+1},beta*N[x]*(int[U[y-x]*I[y],y])}|
\end{equation}


\subsection{Movement}
Movement has to be specified as a pair of transitions,
the local emigration of individuals (which in the simplest case
depends only on local density, but which could be dependent
on neighbourhood density of the same or other species) and
the neighbourhood-dependent immigration of the same individuals.
Technically this preserves only statistical balance rather than
detailed balance of the numbers of individuals, but since we
are deriving equations for the statistics of the population
we can get away with fudging this.  (Another slight fudge is
necessary; you have to tell the package explicitly
which kernels represent movement (see \S\ref{sec:varkernlists}),
to avoid having it think that the immigration rate is
a reproduction term which would require a singular term
in the covariance equations.)

Mathematically:
\begin{equation}
N(x) \to N(x)-1 \mbox{ at rate } -m N(x)
\end{equation}
and
\begin{equation}
N(y) \to N(y)+1 \mbox{ at rate } m \int_\Omega D_m(|y-x|) N(x)
\end{equation}
or in \mma:
\begin{verbatim}
rule5={{{N[x]},{-1},-m*N[x]},{{N[y]},{+1},m*int[Dm[y-x]*N[x],x]}}
MoveKernList={Dm}
\end{verbatim}

\subsection{Exogenous heterogeneity}
In order to put (static) exogenous heterogeneity into the
model, one just has to define the exogenous state variable
as a static variable, and then proceed to use it in the
definitions of other state transitions.  The only trick is
that one has to let \ema\ know that there is no
additional Bernoulli/Poisson variance at small scales (so far this
small-scale variance is required to be zero, this may change in the
future).

For example, a model with density-independent but
spatially heterogeneous mortality rates:
\begin{equation}
\begin{split}
\frac{d\mu(x)}{dt} & = 0 \\
\lim_{r \to 0} \int_0^r c_{\mu \mu}(r') \, dr' & = 0 \\
N(x) \to N(x)-1 & \mbox{ at rate } \mu(x) N(x)
\end{split}
\end{equation}
or
\begin{verbatim}
rule6={{{mu[x]},{0}, 0},{{n[x]},{-1}, mu[x]*n[x]}}
ContVarList={mu}
\end{verbatim}
A straightforward extension in this case
is to have organisms sense the environment, not
just at a particular point, but in a neighborhood
around themselves:
$$
N(x) \to N(x)-1 \mbox{ at rate } N(x) \int H(|y-x|)
\mu(y) \, dy,
$$
or
\begin{verbatim}
|{{n[x]},{-1}, n[x]*int[H[y-x]*mu[y],y]}|
\end{verbatim}
(this model can also be represented by simply convolving
the habitat heterogeneity with the habitat kernel $H$ to
begin with and then treating the result as a static distribution
that is sensed at a particular point \ldots)

I've never built a model with dynamic heterogeneity, but the way to go
about it would be to define a dynamic process that ended up with the
appropriate mean and covariance at equilibrium, then to incorporate it
as a dynamic (rather than a static) state variable.  (This is
essentially what you get with completely asymmetric competition; to
the weaker competitor, the static competitor represents a form of
exogenous spatiotemporal variation.)

\section{Technical details}
\subsection{Variable names}

My conventions are to specify state variables as
single uppercase letters; kernels
as an uppercase letter combined with lowercase letters
or numbers (\verb|Dk| for dispersal, \verb|Mk|
for movement, \verb|Uc|, \verb|Ui|, \verb|U11| etc. for different
competition/infection interactions); spatial variables
as \verb|w|, \verb|x|, \verb|y|, \verb|z|, etc..
I have spelled out Greek letters rather than using
\mma's fancy built-in system for special characters just because
I'm old fashioned.

\textbf{Please keep the following restrictions in mind
when picking variable names --- it's easy to get yourself
in trouble here \ldots}:
\begin{enumerate}
\item{Variables must be legal \mma\ symbols
(e.g. punctuation within names is not allowed);}
\item{Avoid \mma\ reserved words
(\verb|D|, \verb|E|, \verb|I|, \verb|N|, and \verb|O| are
all functions or mathematical constants in \mma);}
\item{Avoid other reserved words that I have
used in this package: \verb|int|, \verb|conv|,
\verb|M|, and \verb|c| should probably
be avoided \ldots (others??)}
\item{\textbf{Position variables}: 
unfortunately, there are fairly strict
conventions on position variables.
I usually take \verb|x| to be the focal site,
\verb|y| to be a site that interacts with the focal site (via
competition, dispersal, etc.), and \verb|z| to be the location
for which covariances are calculated: \verb|v| and \verb|w| are
auxiliary interaction locations (for instance in the case of
density-dependent fecundity where a parent at \verb|y| dispersing its
offspring to \verb|x| may experience competition with site \verb|w| [!])
These conventions are mostly hard-coded, although I have started to
try to relax them. For example,
\verb|x| is automatically set to zero when taking
integrals, to incorporate the effect of spatial homogeneity;
expectations must be taken with the same position variable as
is used in the state transition definition,
or the expected rates will come out zero
(e.g. in the \verb|posvars| argument to \verb|meaneq|
or the \verb|covvars| argument to \verb|coveq|).}
\item{If you plan to translate your equations to \xtc\ 
format (see below), note that \xtc\ is \textbf{not case-sensitive}
(so \verb|f| and \verb|F| will get confused), and may have
some of its own reserved words (\verb|exp|, \verb|gamma|, etc.).}
\end{enumerate}

\subsection{Miscellaneous}

\verb|TableForm[model]| will give a prettier display of the
model definitions.  (There are some sketchy functions to translate
state variables etc. to \LaTeX\ format.)

\section{Example so far}
Now that you've gotten this far, here's a very simple example:
the rules for the ``spatial logistic equation'' with density-dependent
mortality.
In this model, we assume that individuals experience density-independent
mortality at \emph{per capita} rate $\mu$ (i.e., the probability of
an individual dying in a small interval $\Delta t$ is $\mu \Delta t$);
they experience an additional density-dependent mortality,
linearly increasing with slope $\alpha$ as 
a function of local density weighted by a competition
kernel $U$; and they reproduce at constant $\emph{per capita}$
rate $f$, after which their offspring (seeds, whatever) immediately
disperse away from the parent with dispersal kernel $D$.

In mathematical form:\\
\begin{tabular}{lll}
{\bf transition} & {\bf rate} & {\bf description} \\
$N(x) \to N(x)-1$ & $\mu N(x)$ & density-independent mortality \\
$N(x) \to N(x)-1$ & $\alpha N(x) \int U(y-x) N(y) \, dy$ & density-dependent mortality \\
$N(x) \to N(x)+1$ & $f \int D(y-x) N(y) \, dy$ & reproduction and local dispersal
\end{tabular}

Enter the model in \mma\ like this:
\begin{verbatim}
spatlogist = {{{N1[x]},{-1},mu*N1[x]},
              {{N1[x]},{-1},N1[x]*alpha*int[U[y-x]*N1[y],y]},
              {{N1[x]},{+1},f*int[D1[y-x]*N1[y],y]}}

\end{verbatim}

And here's what a little piece of the \mma\ session looks like:
\begin{verbatim}
In[5]:= In[5]:= TableForm[spatlogist]

Out[5]//TableForm= N1[x]   -1   mu N1[x]

                   N1[x]   -1   alpha int[N1[y] U[-x + y], y] N1[x]

                   N1[x]   1    f int[D1[-x + y] N1[y], y]

\end{verbatim}

A few notes:
\begin{itemize}
\item{you can combine density-dependent and density-independent
mortality in a single rule, but it seems clearer if you separate them;}
\item{avoid the temptation to use {\tt N} as a state variable or
{\tt D} as a dispersal kernel ---
they're both reserved words in \mma.}
\item{\mma\ sorts variables in expressions alphabetically (so, 
for example, the kernel {\tt D1} comes before the state variable
{\tt N1} in the integral, but {\tt U} comes after).}
\end{itemize}

\section{Generating moment equations}

\subsection{Variable and kernel lists}
\label{sec:varkernlists}
The first step is to extract lists of state variables, spatial
position variables, and kernels from the
equations (\emph{this should now be done completely automatically}).
State variables are extracted from the first entries in each transition definition;
kernels are identied as any \mma\ function in the transition
rates that is neither an integral nor a state variable.)
If you should need to do this manually for some reason,
\begin{equation}
\verb|VarList=getvarlist[model]|,
\end{equation}
\begin{equation}
\verb|KernList=getkernlist[model]|,
\end{equation}
and
\begin{equation}
\verb|XVarList=getxvarlist[model]|
\end{equation}
will extract the variables, kernels, and position variables respectively.

In order to do the bookkeeping correctly, movement kernels do
need to be specified separately by setting 
e.g. \verb|MoveKernList={Mk}|, and
continuous (non-discrete) environmental
variables need to be specified as
part of \verb|ContVarList|.

\subsection{Mean and covariance equations}
Once you've defined your model you can
generate the moment equations simply by saying
\begin{equation}
\verb|meaneqs[model]|
\end{equation}
to get the equations for time-evolution
of the mean densities (RHS of ordinary differential
equations for each state variable's density in order of
\verb|VarList|, which is alphabetical).
Mean densities of state variables are denoted
by \verb|M[n]|.  The expression \verb|cbar[U,m,n]|
denotes an kernel-weighted average covariance,
corresponding to $\int U(r) c_{mn}(r) \, dr$.
``Naked'' covariances $c_{mn}$ should not appear in
the mean equations!

Use \begin{equation}
\verb|coveqs[model]|
\end{equation}
for the covariance equations (RHS of integral equations
in pairwise order).  The spatial lag variable $r$ is
implicit throughout these equations.
Covariances of state variables are \verb|c[m,n]|;
convolutions between covariances and kernels
are \verb|conv[c[m,n], K, z]| (= $(K*c)(r) = \int U(s) c_{mn}(r+s) \,
ds$). (FIXME:
The spatial variable \verb|z| is redundant/implicit and
should be dropped. \verb|implicitvar|?)
Third central moments (which always appear inside an integral,
\verb|int|) \verb|M3| depend on three state variables.

If you want to see exactly what \mma\ is doing to reduce
your model to moment equations, try setting
\verb|exdebug=True| and running \verb|meaneqs| or \verb|coveqs|
again \ldots \verb|intdebug=True| should provide information
on the integration rules.

\subsection{Moment closure}

\fixme{closeM3 is now closeM!}

Third central moments \verb|M3| have to be dealt with before
you can proceed to analyze or numerically integrate the equations.
The function for performing moment closure on a set of equations
is \verb|closeM3|.  It has a series of options that can be
set; the most important is \verb|closure|.
\begin{itemize}
\item{To set third central moments to zero (power-1 closure),
use \verb|closure->power1|.}
\item{To use an asymmetric power-2 closure (where we decompose
the conditional probability of a triple $XYZ$ into $p_{XYZ} =
p_{XY} p_{YZ}/p_Y$), use
\verb|closure->power2a|.
There is an additional option for \verb|power2a|,
\verb|closepos| (default \verb|0|), which specifies the spatial location
of the ``middle'' of the triple.}
\item{symmetric power-2 closure uses
$$
p_{XYZ} = \frac{1}{2} \left(
\frac{p_{XY} p_{YZ}}{p_Y} +
\frac{p_{XZ} p_{YZ}}{p_Z} +
\frac{p_{XY} p_{XZ}}{p_X} -
p_{X} p_{Y} p_{Z} \right).
$$
The last term inside the parentheses is necessary for
consistency (REF Dieckmann and Law, IIASA book).}
\item{weighted power-2 closure is a generalization of the preceding
two closure, due to Murrell, Law and Dieckmann:
$$
p_{XYZ} = \frac{1}{\alpha+\gamma} \left(
\alpha \frac{p_{XY} p_{YZ}}{p_Y} +
\beta \frac{p_{XZ} p_{YZ}}{p_Z} +
\gamma \frac{p_{XY} p_{XZ}}{p_X} -
\beta p_{X} p_{Y} p_{Z} \right).
$$
There is an extra parameter, {\tt closeweights}, which
defaults to $\{\alpha,\beta,\gamma\}=(4,1,1)$: 
{\tt closeweights}=(1,0,0) corresponds to {\tt power2a} and
{\tt closeweights}=(1,1,1) corresponds to {\tt power2s}.}
\item{To do power-3 closure, making the assumption that
$p_{XYZ} = \frac{p_{XY} p_{YZ} p_{ZX}}{p_X p_Y p_Z}$:
use \verb|closure->power3|.}
\end{itemize}

\verb|closeM3| also handles what I call ``third-moment singular
terms'', which arise when the triple $XYZ$ is really just a pair
($X=Z$, etc.) or a single location ($X=Y=Z$).  These third-moment
singular terms are taken into account if \verb|singm->True|, which is
the default; set \verb|singm->False| to ignore them. (Since a typical
third-moment singular term $S$ appears in the form $K(r) S(r)$ rather
than the form $\int K(r) S(r) \, dr$, one can argue that
third-moment singular terms are of the same order as the scale
parameter and can be neglected when all kernels are relatively
long-range: Bolker and Pacala 1999 make this assumption.)

Singular terms can only arise between a pair $X(x) Y(y)$ in a triplet
$X(x) Y(y) Z(z)$ if (1) the state variables $X$ and $Y$ are identical
(otherwise the locations $x$ and $y$ could not refer to the same
point, because it would have to be occupied by two different
individuals) and (2) $x$ and $y$ are not the lag variables involved in
the current covariance calculation (otherwise this singular term would
only appear in the covariance at lag zero, which we do not calculate
anyway).  The option \verb|covpos| (default \verb|{0,z}|) specifies
the lag variables for the current covariance calculation, so that
criterion 2 can be determined.

It is often useful to simplify the equations further
after applying moment closure rules: something like
\begin{verbatim}
simpeqs = eqs //. intrules
\end{verbatim}
will apply a set of rules for simplifying integrals.
There is an important ``state variable'' called
\verb|singtaken|; if \verb|True|, it tells
\mma\ that singular terms have already been taken out of
all covariances and third moments (so it doesn't try to extract them
again).
\verb|closeM3| has a \verb|setsingtaken| option (\verb|True|
by default) that sets \verb|singtaken=True|.
Set \verb|singtaken| back to \verb|False| if you want to derive
more moment equations from transition rules in the same \mma\ session.
\fixme{perhaps \code{meaneqs}/\code{coveqs} should automatically
toggle \code{singtaken} back to \code{False}?}

At this point you may also want to use the function
\verb|implicitvar|, which drops a lag variable 
(by default \verb|z|) from the equations
and makes it implicit.  For example, \verb|c[s,s,z]| (which
may already have been implicitly shortened from \verb|c[s,s,0,z]|) gets
shortened to \verb|c[s,s]|.  Once you reach this stage there should be
no explicit spatial position variables left in the equation.

\subsection{Scaled-covariance and correlation equations}
The equations produced at this point are equations for the
covariance density.  You may prefer to work with the equations
in another scaling, particularly if you plan to analyze invasion
criteria.

Here are some of the possible scalings one could use (not all of which
are available in this package): \\
\begin{tabular}{ll}
probability & $p_{xy}$ \\
conditional probability & $p_{x|y} = \frac{p_{xy}}{p_y}$ \\
covariance & $C_{xy} = p_{xy} - p_x p_y$ \\
scaled covariance & $ SC_{xy} = \frac{C_{xy}}{p_x} = \frac{p_{xy}-p_x p_y}{p_x} = p_{y|x} - p_y$ \\
correlation & $ c_{xy} = \frac{C_{xy}}{p_x p_y} = \frac{p_{xy}-p_x p_y}{p_x p_y} = \frac{p_{xy}}{p_x p_y} - 1$ \\
``multiplicative covariance'' (?) & $\frac{p_{xy}}{p_x p_y} = c_{xy}+1$
\end{tabular} 

\emph{8 June 2002: changed order/definition of scaled covariance so that $SC_{xy}$ 
is now proportional
to neighborhood density/conditional probability of $y$ near $x$, for consistency
with competition coefficients, so that $\alpha_{xy} \overline{SC_{xy}}$ is the
competitive effect of $y$ on $x$ \ldots}

\begin{verbatim}
covtocor[coveqs,Teqs]
\end{verbatim}
translates a set of covariance equations \verb|coveqs|
derived from an original set of transition probabilities
\verb|Teqs| to correlation equations.  (You need to specify
the original transition equations again because the
translation involves the expressions for changes in the
mean densities as well as covariances.)

Similarly,
\begin{verbatim}
covtoscov[coveqs,Teqs]
\end{verbatim}
translates to scaled covariances.  Note that there will
be more scaled-covariance equations than covariance equations
($n^2$ rather than $n(n+1)$),
because scaled covariances are not symmetric
($SC_{xy}=p_x/p_y SC_{yx}$).
You will also need something like
\begin{verbatim}
meqs /. covcorrules
\end{verbatim}
to translate any covariance (or averaged covariance)
terms in the mean equations to scaled equivalents.

\section{Analyzing models}
You can take your equations at this point and do whatever you like
with them: analyze equilibria, quasi-equilibria, integrate the
equations numerically, etc..
If you choose to continue the analysis within \mma,
there are various \mma\ commands (built-in or included in the \ema) that may come in handy.
These are just a few hints: for complete details on analysis
you'll have to read the various \mma\ documentation.
\begin{itemize}
\item{Use the pattern rule \verb|/.| to
substituting simplifying combinations of parameters or
numerical values of parameters.
To apply multiple rules you may want to use \verb|Flatten|.
Examples:
\begin{verbatim}
simpeqs = Simplify[eqs /. {h->(rho-d-r), alpha->A-m}];
simpeqs = eqs /. {a -> 1.3, b-> 2.7}
\end{verbatim}}
\item{\verb|fourtrans|: apply a symbolic Fourier transform
to covariance equations (all this function
does is to transform convolutions into products and
\emph{vice versa}; it's only really worth it if
the only covariance terms in your equations are in terms
of $\bar c$ terms or convolutions, in which case you will
only have products in the Fourier-transformed equations).}
\item{\verb|zerocbar|: drop spatial terms in mean equations,
leading to the mean field model}
\item{\verb|sameK|: set all kernels equal to each other}
\item{Substitute functional forms of kernels (in 
space or frequency domain), e.g. \verb|/. K -> 1/(1+q^2)|}
\item{substitute mean-field solutions in covariance equations
(approximate solution when effects of space are small $\to$
densities are close to MF values)}
\item{\verb|momcollect|: collect covariance (and possibly
third-moment) terms in a set of equations.}
\end{itemize}

\subsection{Equilibrium/quasiequilibrium solutions}
Integration of {\tt cbar}, etc..

\subsection{Kernel identities}
If we assume Laplacian kernels ($K(r) = \lambda/2 e^{-\lambda|r|}$)
in one dimension or modified-Bessel kernels ($K(r) = K_0(\lambda r)$)
in two dimensions, we can get some mileage out of a set
of \emph{kernel identities}.  What is special about these
kernels is that they have the same convenient functional
form in the frequency domain, $\lambda^2/(\lambda^2+q^2)$.
``Rational functions'' of kernels (e.g., of the form $K_1/(c+K_2)$)
in the frequency domain can be reduced by partial fractions to
somewhat simpler forms, and I've encoded many of these rules
in \mma (see {\tt kern1.ma} [which probably needs significant cleanup
and documentation before anyone else could bear to look at it]).

\subsection{Translation to \xtc}
The package can output any equations you have derived as an input file
for \xtc, a program by Bard Ermentrout that does numerical
integration of one-dimensional integral equations
(URL \url{http://www.math.pitt.edu/~bard/bardware/xtc-1.3.tar.Z}).  
Many of the
combinations of moment closure/Fourier transform, along with suitable
assumptions about radial symmetry, allow the moment equations to be
written as 1D integral equations.  Of course, any equations for an
explicitly 1-D system can be written as 1-D equations.  In addition,
the Fourier transformed version of MEs with power-1 closure or the
regular (space domain) version of MEs with power-2A closure come out
as 1-D equations, with the added bonus that the equations contain no
explicit convolutions, which are possible, but slow, in \xtc\ 
(although you may have to drop third-moment singular terms to make this
work).

I had written a whole set of C codes for numerically integrating
moment equations, but I have turned to \xtc\ because it is more
flexible.  I have fixed some bugs in \xtc\ and added some 
functionality (batch mode, ability to read in starting conditions
for continua from a file); for the time being, you can get my
version (latest version is 1.3c) from my web page.
In order to run \xtc\ on Windows,
you will have to compile it (it comes as fairly portable C code) and
get yourself an X-windows emulator to run the graphics (see \ldots).
I intend to write a version of \xtc\ that can be compiled with graphics
turned off (and run only in batch mode), but this hasn't happened yet.
Ermentrout's {\tt xppaut} package, for numerical analysis of ODEs, is
far slicker and has had a lot of work in the past few years.  He's
adapted it to solve most of the PDE-like problems he is interested in,
so upgrades for \xtc\ (which is written and distributed entirely at
his own initiative and expense) are less likely to happen.
Nevertheless, I find \xtc\ sufficiently powerful and convenient that I
prefer it to {\tt xppaut} (or to my own code) for these applications.

The command \verb|mmatoxtc[meqs,coreqs,parvals,cparvals,mstart,opts,xtcopts,backup]|
produces an \xtc\ input file with a given set of mean equations,
covariance/scaled covariance/correlation equations, parameter values,
continuum parameters,
starting values for the mean densities, and options.
Parameters are expressed in the form \verb|{p1->1.0,p2->2.0}|, etc.;
kernel parameters have to be expressed as \verb|KP1| for the first
parameter for kernel \verb|K|, etc..
Remember that \xtc\ is not case sensitive, so choose your variable
names with care.
Options are expressed as \verb|option->value|: your choices are
\begin{itemize}
\item{dimension, \verb|dim| (1 or 2);}
\item{Fourier-transformed equations, \verb|ftrans| (\verb|True| or
\verb|False|);}
\item{whether the equations are for scaled covariances rather
than covariances (in which case $c_{12}$ and $c_{21}$ are not
identical): \verb|sccov=False| by default}
\item{output file, \verb|outfn| (a quoted file name: note that
the code automatically makes numbered backups if you try to
overwrite an existing file);}
\item{total arena length (\verb|length=10| by default)
and grid scale (\verb|griddx|=0.02 by default);}
\item{kernel shapes:
defaults are 
\begin{itemize}
\item{\verb|cauchy|=$\lambda^2/(\lambda^2+q^2)$ if
Fourier-transformed}
\item{\verb|mexp|=$\lambda/2 \exp(-\lambda x)$ (back-to-back exponential, renamed to avoid conflicting
with \xtc's own exponential function) in 1 dimension}
\item{or \verb|mexp2d|=$\lambda^2/(2 \pi) \exp(-\lambda r)$ in 2
dimensions}.
\end{itemize}
The other choices are
\begin{itemize}
\item{\verb|gauss|($\sigma$)=$1/(\sqrt{2 \pi} \sigma) \exp(-(x/\sigma)^2/2)$}
\item{\verb|gauss2d|($\sigma$)=$1/(2 \pi \sigma^2) \exp(-(r/\sigma)^2/2)$}
\item{\verb|beta|($a$,$b$)=$H(1-x) \frac{\Gamma(a+b)}{\Gamma(a)\Gamma(b)} x^{a-1} (1-x)^{b-1}$}
\item{\verb|sbeta|($a$,$b$,$s$) [scaled beta]=$H(s-x) \frac{\Gamma(a+b)}{\Gamma(a)\Gamma(b)} (x/s)^{a-1} (1-(x/s))^{b-1}$}
\item{\verb|errdist|($\tau$,$d$)=$\frac{\Gamma(2/\tau)}{2\tau d (\Gamma(1+1/\tau))^2}
\exp(-(x \Gamma(2/\tau)/(d \tau \Gamma(1+1/\tau)))^\tau)$ (a
``generalized Gaussian'' distribution with mean dispersal distance $d$
and exponent $\tau$: proportional to exponential with $\tau=1$,
Gaussian with $\tau=2$)}
\end{itemize}
}
\item{starting values of covariances, \verb|cstart|}
\item{various \xtc\ settings (as a list of \verb|xtcopts|):
\verb|deltat|, \verb|tfinal|, \verb|nout|,
\verb|method|, \verb|fast|, \verb|bound|, \verb|bufsize|,
\verb|transient|;
you can also set these within \xtc\ (see \xtc\ documentation for details).}
\end{itemize}

For example:
\begin{verbatim}
mmatoxtc[meqs,ceqs2,{f->5,mu->1,alpha->1,UP1->1,DkP1->1},{},{0.1},0,
{dim->1,ftrans->False,outfn->"DDM1a.xtc"},{nout->500}]
\end{verbatim}

You will have to read the documentation of \xtc\ for details
of how to run and analyze the numerical models,
but the basic steps are:
\begin{itemize}
\item{\verb|xtc foo.xtc| at a shell prompt;}
\item{{\bf G} for Go, to run the equations;}
\item{{\bf W}indow/{\bf F}it 3d to adjust colors in the plot;}
\item{{\bf 2}D graphics/{\bf 4} (U vs. t) to bring up a dialog
box for displaying scalar quantities (densities or $\bar c$ values),
and {\bf W}indow/f{\bf I}t 2d to adjust/display the plot;}
\item{{\bf R}edraw menu to redraw 2- or 3d plots after they are
changed or obscured by other windows;}
\item{{\bf F}ile/{\bf S}ave 3d to save data shown in the 3d plot,
or {\bf F}ile/s{\bf A}ve 2d for the 2d plot;}
\item{{\bf F}ile/{\bf Q}uit/{\bf Y}es to stop.}
\end{itemize}

You can also call sets of batch runs from \xtc.
\begin{itemize}
\item{\verb|runxtc| is the basic function for running moment
equations in \xtc\ from \mma.  It writes out an output file,
calls the operating system to run \xtc, and returns the results
(\emph{describe format}).}
\item{\verb|findxtceq0| runs equations to equilibrium,
returning the last two time steps}
\item{\verb|findxtceq| runs equations to equilibrium,
checks for convergence, and returns the state at the last
time step} 
\item{\verb|findxtcinv| runs an invasion analysis in \xtc,
running one species to monoculture equilibrium and then
introducing the other at low densities.
\textbf{Parameters:} 
\begin{itemize}
\item{mean equations}
\item{covariance equations}
\item{parameter values}
\item{starting values for mean (list of same length as mean equations)}
\item{starting values for covariances}
\end{itemize}}
\item{\verb|findxtceqtab| finds the equilibrium for a series
of parameter values (specify the parameters as a list of sets
of parameters)}
\item{\verb|findxtceqstab| finds the equilibrium for a series
of starting values}
\end{itemize}
The order of variables in the output is: 
(time, mean 1, \ldots, mean n, \{ cov 1 \}, \ldots, \{ cov n \},
cbar 1, \ldots, cbar n). 
Note that \verb|findxtceqtab| goes to some trouble to start each
new run from the \emph{ending} condition of the previous run, so
that equilibrium should be achieved faster (if changes in parameters
from run to run are fairly gradual).

Some other notes about running moment equations under \xtc:
\begin{itemize}
\item{the Fourier-mode and 2D calculations are not tested yet}
\item{length is set in two places in the \xtc\ file (as an
option and as a parameter), be careful to change them both
if you are editing the file}
\item{You do have to be careful when numerically integrating moment
equations:
\begin{itemize}
\item{power-1 closures will often blow up if you start the mean too
small (this is the fault of the closure itself, not the numerical
methods)}
\item{if you don't set the length of the arena large enough, you may
have problems because of the way \xtc\ extrapolates beyond the edge of
the arena.  If the value of the covariance at the maximum distance is
not close to zero, watch out}
\end{itemize}}
% \item{at the moment it isn't particularly easy to simulate
% invasions to a monoculture equilibrium under \xtc.  I use the {\tt
% heav} built-in Heaviside function ($H(x)=0$ if $x\leq 0$, otherwise 
% $H(x)=1$) to simulate a transient: the derivatives of the invader,
% the invader autocovariance, and the cross-covariance are multiplied by
% {\tt heav(t-trans)} so that they can't increase or decrease until
% $t>\mbox{trans}$.  This is not exactly right since the resident still
% has to deal with a small, fixed, randomly dispersed
%($c_{II}=c_{IR}=0$) population of the invader even before the invasion
%occurs, but it works pretty well}
\item{you can't use an expression within a convolution --- this is OK
except for multiple convolutions (e.g. $(Uc)*c$ or $(U*c)*c$), can
be solved by adding another continuum fixed variable that defines 
the expression}
\item{\xtc\ wish list:
\begin{itemize}
\item{zoom/thinning capability for display of 3D surfaces}
\item{mixed extensions (e.g. even on left boundary, zero on right
boundary)}
\item{label 2-D graph window with file}
\item{direct output (png?) of 3D color graph}
\item{more robust generally (e.g. segfaults if a scalar
equation is undefined; crazy if one tries to use expressions within
a convolution); better error reporting?}
\item{ability to define a default 2D graph}
\end{itemize}}
\end{itemize}

\section{Wish list/notes}
\begin{itemize}
\item{what changes would actually be required for non-infinitesimal patches?
Could you still assume Poisson statistics at a patch or would you have
to make zero-dimensional (distributional) moment closure and track
variance separately?}
\item{need some rules for environmental variables with non-Poisson variances:
this does come up occasionally (e.g. when organisms sense the environment
through a ``habitat kernel'')}
\item{directional dispersal??}
\item{Note order of covariances etc.: alphabetical}
\item{Further functions for quasi-eq/invasion analysis??
(some of these already exist)}
\item{Fix/sort out sensitivity to spatial positions}
\item{correlated dispersal??}
\item{correlated disturbance??}
\item{tests: power-1, power-2, power-3 closure, 1- and 2-D,
with and without singular terms, space and frequency domain, DDM/DDE/DDF.
Cf. simulations \ldots different parameter regimes \ldots}
\item{dealing with self-competition: \verb|U[0]->0|;
works OK but is wrong in \xtc\ output if we change to FT
kernels and leave $U(0)$ unchanged}
\item{FIXME: order of n1, n2 in Cbar variables is still getting switched around
by mmatoxtc}
\item{FIXME: change variables when taking Fourier Transform?}
\item{FIXME: should \verb|closeM3| automatically apply \verb|intrules|?}
\item{tests:
\begin{itemize}
\item{replicate equations/numerical solutions in Bolker \& Pacala
papers;}
\item{replicate equations/numerical solutions in Bolker '99;}
\item{replicate equations in Brown \& Bolker ms.;}
\end{itemize}}
\item{how many possible models/numerical sims are there?
3 closures (power-1, power-2A, power-3) $\times$ 
singular terms (1,2) $\times$ dimension (1,2) $\times$
state variables (covariances, scaled cov., correlations)
$\times$ number of species (1,2) $\times$ competition type
(DDM, DDE, DDF) $\times$ special cases (neutral, symmetric,
equal-scales, quasi-equilibrium, etc. etc.) \ldots}
\end{itemize}

\section{Examples}
\subsection{Disturbance}
\section{Tests}
\section{History}
\subsection{Changes to code:}

\subsection{Changes to documentation:}
\begin{itemize}
\item{July 2004: tweaking.  Eliminated
``tests'' section (until I can get it a bit cleaned up)}
\end{itemize}

\section{Notes}
6 July 2004: 
\begin{itemize}
\item{had some trouble with order of variables
in an M3 term (not recognized) --- had to use
{\tt resort}.  Put in automatically at some point?}
\item{must use \verb|resort[c1,3] //.intrules| before
trying to use \verb|closeM|}
\end{itemize}


\end{document}


\begin{itemize}
\item{Model:
\begin{itemize}
\item{spatial logistic: 
\begin{itemize}
\item{DDM, DDE, DDF}
\item{self-competition?}
\item{clustered dispersal?}
\item{1-sp., 2-sp.}
\end{itemize}}
\item{pred-prey/epidemic:
\begin{itemize}
\item{self-limitation in pred/prey?}
\item{movement by pred/prey?}
\item{efficiency: $<1$, $=1$, $>1$}
\end{itemize}}
\end{itemize}}
\item{Derivation:
\begin{itemize}
\item{Closure: power-1, 2a, 2, 3}
\item{Singular terms?}
\item{Spatial state variable: covariance, scaled cov., correlation, cond. prob.?}
\end{itemize}}
\item{Arena:
\begin{itemize}
\item{Dimension (1,2)}
\item{Exogenous heterogeneity (in mortality, fecundity, movement, predation,
competition)?}
\end{itemize}}
\item{Numerical solution:
\begin{itemize}
\item{Real (1D)}
\item{Fourier (1D, 2D)}
\end{itemize}}
\end{itemize}

Default \xtc\ params. (equilibrium):
\begin{verbatim}
delta_t=0.025, tfinal=0.1, nout=500, method=0, fast=1, bound=1000,
length=10, grid=500
\end{verbatim}

Default sim params (equilibrium):
\begin{verbatim}
trans=0.0, end_t=50.0, tau=0.02, t_rpt=0.02,
x_rpt=0.1, ranseed=1001, nplantbox=1, nsumbox=51,
len=50, covdist=50.0, pdedist=50.0 
\end{verbatim}

\begin{tabular}{ccccp{2in}c}
\textbf{Model} & \textbf{Method} & \textbf{Equations} & \textbf{Parameters}
& \textbf{Notes} & \textbf{Results} \\
spatial logistic (DDM) &
power-1/1D/real & & $f=5$, $\mu=1$, $\alpha=1$, $\lambda_U=\lambda_D=1$ &
no self-comp; start=(0.1,0)
[start=0.01 crashes PDE;
$f=2$ makes nbrhd size too small]
(eqns check with BP TPB 1997)
power-1 more sensitive with scaled cov: need
higher starting cond. (n=0.4)
Gear not necessary despite initial red herrings 
\end{tabular}
\begin{verbatim}
DDM1a (d1=u1=1)
SIMULATION:
           mean        cbar      nbrhd
mvec  3.9105953 0.157937418 3.95175036
sdvec 0.0467677 0.080014386 0.04818047
lovec 3.8189306 0.001109221 3.85731663
hivec 4.0022600 0.314765614 4.04618409

power-1:    3.944326 0.219598 4.000000
power-1 (grid=100):
            3.944117 0.220408  
(R eq/100:) 3.937929 0.244425 3.9999985
power-2a:   3.553164 1.587681 4.000000
power-2s:   3.566640 1.545638 4.000000 
power-3:    3.949307 0.200201 4.000000

power-1/sc: 3.944047 0.055953 4.000000 
power-2a/sc: 3.553164 0.446836 4.000000 
power-2s/sc: 3.583591 0.416409 4.000000
power-3/sc:  3.949307 0.050693 4.000000  

(d1=u1=2)
            mean      cbar      nbrhd
mvec  3.84410204 0.3634281 3.93921978
sdvec 0.05841351 0.1140862 0.07204214
lovec 3.72961156 0.1398191 3.79801719
hivec 3.95859252 0.5870370 4.08042237

power-1:    3.888226 0.434604 4.000000 
power-2a:   3.190993 2.581536 4.000000 
power-2s:   3.393463 2.059003 4.000219 
power-3:    3.907608 0.361032 4.000000 

power-1/sc:  3.887425 0.112575 4.000000 
power-2a/sc: 3.190993 0.809007 4.000000 
power-2s/sc: 3.28206  0.71794  4.00000 
power-3/sc:  3.906392 0.093608 4.000000 
\end{verbatim}

Quasi-equilibrium scaled average covariance ($n=0$):
parameters as above but $\Delta t=0.001$,
reporting time 0.10, 500 steps ($t=5$).
\begin{verbatim}
power-1: blows up
power-2a: 0.983557
power-2s: 1.107480
power-3: infinite derivative
\end{verbatim}
Sim results for quasi-equilibria??

Results for heterogeneous landscapes:
($f=5$, $\bar\mu=1$, $\alpha=1$, $\lambda_U=\lambda_D=1$,
mortality landscape is Gaussian with $\sigma=1$: \\
\verb|transient=5, nout=1, grid=250, delta_t=0.025, tfinal=0.1|) \\
\begin{tabular}{lcccc}
closure  & $\bar n$ & $\bar c_{nn}$ & $c_{\mu n}(0)$ \\
power-1  & 3.98446 & 0.404548 & -0.342624 \\ 
power-2a & 3.55322 & 1.5875 & 1.66533 $\times 10^{-17}$ \\
power-2s & 3.5224  & 1.79914 & -0.116889  \\
power-3  & 3.98923 & 0.340262 & -0.297308
\end{tabular}

Calculate nbrhd size?
(for spatial logistic: R (f/mu), nbrhd size (K/d1), d1/u1 completely determine model
[r=time scale, K=spatial scale])

\section{Examples}
\section{History}
\subsection{Changes to code:}
\begin{itemize}
\item{July 2004: tweaked to fix addition of {\tt debug}
flag
\item Jan 2009: added \code{cq} to momcollect
\end{itemize}

\subsection{Changes to documentation:}
\end{document}
