camodel {chouca}

Definition of a stochastic cellular automaton

Description

High-level definition of a stochastic cellular automaton

Usage

camodel(
  ...,
  neighbors,
  wrap,
  parms = list(),
  all_states = NULL,
  check_model = "quick",
  verbose = FALSE,
  epsilon = sqrt(.Machine[["double.eps"]]),
  fixed_neighborhood = FALSE
)

transition(from, to, prob)

Arguments

...

A number of transition descriptions, as built by the transition function (see Details and Examples)

neighbors

The number of neighbors to use in the cellular automaton (4 for 4-way or von-Neumann neigborhood, or 8 for an 8-way or Moore neighborhood)

wrap

If TRUE, then the 2D grid on which the model is run wraps around at the edges (the top/leftmost cells will be considered neighbors of the bottom/rightmost cells)

parms

A named list of parameters, which should be all numeric, single values

all_states

The complete set of states of the model (a character vector). If unspecified, it will be guessed from the transition rules, but it is a good idea to pass it here to make sure the model definition is correct.

check_model

A check of the model definition is done to make sure there are no issues with it (e.g. probabilities outside the [0,1] interval, or an unsupported model definition). A quick check that should catch most problems is performed if check_model is "quick", an extensive check that tests all possible neighborhood configurations is done with "full", and no check is performed with "none".

verbose

Whether information should be printed when parsing the model definition.

epsilon

A small value under which the internal model coefficients values are considered to be equal to zero. The default value should work well here, except if you run models that have extremely small transition probabilities (<1e-8).

fixed_neighborhood

When not using wrapping around the edges (wrap = FALSE), the number of neighbors per cell is variable, which can slow down the simulation. Set this option to TRUE to consider that the number of neighbors is always four or eight, regardless of the position of the cell in the landscape, at the cost of approximate dynamics at the edges of the landscape.

from

The state from which the transition is defined

to

The state to which the transition is defined

prob

a one-sided formula describing the probability of transition between the two states (see Details section for more information).

Details

This help page describes in detail technical points related to the definition of models in chouca. If this is your first time working with chouca, you may like the longer introduction in the vignette, accessible using vignette("chouca-package").

camodel allows defining a stochastic cellular automaton model by its set of transition rules. These are defined by a set of calls to the transition() function. Each of these calls defines the two states of the transition, and the probability as a one-sided formula involving constants and the special vectors p and q.

transition() calls takes three arguments: the state from which the transition is defined, the state to which the transition goes, and a transition probability, defined as a one-sided formula. This formula can include numerical constants, parameters defined in the named list parms, and any combination of p['a'] and q['b'], which respectively represent the proportion of cells in a landscape in state 'a', and the proportion of neighbors of a given cell in state 'b' ('a', and 'b' being, of course, any of the possible states defined in the model). Such formula could typically look like ~ 0.2 + 0.3 * p["a"] + q["b"]. See below for examples of model definitions.

It is important to remember when using this function that chouca only supports models where the probabilities depend on constant parameters, the global proportion of each state in the landscape, and the local proportion of cells around a given cell. In other words, all transition probabilities should have the following functional form:

a_0 + \sum_{k=1}^S g_k(q_k) + s(q, q) + s(p, q) + s(q, q)

where a_0 is a constant, g_k are univariate functions of q_k, the proportions of neighbors of a cell in state k, and q is the vector containing all the q_k for k between 1 and S, the total number of states in the model. Similarly, p is the length-S vector containing the proportion of cells in each state in the whole grid. s above is the sum, defined for two vectors x = (x_1, ..., x_S) and y = (y_1, ..., y_S) as

a_1 x_1^{\alpha_1} y_1^{\beta_1} + a_2 x_1^{\alpha_2} y_2^{\beta_2} + a_3 x_1^{\alpha_3} y_3^{\beta_3} + a_4 x_2^{\alpha_3} y_1^{\beta_3} + a_4 x_2^{\alpha_3} y_2^{\beta_3} + \dots + a_K x_S^{\alpha_K} y_S^{\beta_K}

where the a_k, \alpha_k and \beta_k are constants for all k, and K is the total number of terms (equal to S^2). Note that \alpha_K and \beta_K are capped to 5. This can be overriden using options(chouca.degmax = n), but we do not recommend changing it as higher values typically make the package slow and/or leads to numerical instabilities. The functions g_k above can be any univariate functions of q_k, so chouca effectively supports any type of transition rule involving the neighborhood of a cell, including some 'threshold' rules that involve a single state (and only one). For example, a rule such as "more than 5 neighbors in a given state make a cell switch from state A to B" is OK, but combining states may not be supported, such as "more than 5 neighbors in state A *and* 2 in state B means a cell switches from A to B". When in doubt, just write your model, and chouca will tell you if it cannot run it accurately by running model checks.

Model checks are controlled by the argument check_model. When set to "quick" or "full", a check is performed to make sure the functional form above is able to accurately represent probabilities of transitions in the model, with "full" enabling more extensive testing, and "none" removing it entirely. Coefficients in the formula above are rounded down to zero when below epsilon. This may be an issue if your transition probabilities are close to zero: consider reducing epsilon to a smaller value in this case, or adjusting your model parameters.

When space does not wrap around (wrap = FALSE), cells in the corners or in the edges will have a lower number of neighbors. The proportions of cells in a given state k, q_k, will thus be computed with a reduced number of cells. For example, a cell in a corner will have only 2 neighbors when using a 4x4 neighborhood, so q_k is computed using only two cells, and can be only equal to 0, 0.5 or 1.

To run a model once it is defined, the function run_camodel can be used, or run_meanfield for a mean-field approximation. An initial landscape for a simulation can be created using generate_initmat.

You can update a model definition with new parameters (all of them or a subset) using the update method. The model graph with the different states and transitions can be displayed using the plot method (this requires the package igraph).

Value

This function returns a list object with class ca_model, with the following named components. Please note that most are for internal use and may change with package updates.

transitions

the list of transitions of the model, as returned by transition

nstates

the number of states of the model

parms

the parameter values used for the model

beta_0,beta_q, beta_pp, beta_pq, beta_qq

internal tables used to represent probabilities of transitions when running simulations, these tables are for internal use and probably not interesting for end users, but more information is provided in the package source code

wrap

Whether the model uses a toric space that wraps around the edge

neighbors

The type of neighborhood (4 or 8)

epsilon

The epsilon values used in the model definition, below which transition probabilities are assumed to be zero

xpoints

(for internal use only) The number of values used to represent the proportion of neighbors of a cell in each state

max_error, max_rel_error

vector of numeric values containing the maximum error and maximum relative error on each transition probability

fixed_neighborhood

flag equal to TRUE when cells have a fixed number of neighbors

Functions

• transition():

See Also

run_camodel, run_meanfield, generate_initmat, run_meanfield, update.ca_model, ca_library

Examples

# Redefine Kubo's 1996 forest gap model
kubo <- camodel(
  transition(from = "TREE",
             to   = "EMPTY",
             prob = ~ d + delta * q["EMPTY"] ),
  transition(from = "EMPTY",
             to   = "TREE",
             prob = ~ alpha),
  parms = list(d = 0.125,
               delta = 0.5,
               alpha = 0.2),
  all_states = c("EMPTY", "TREE"),
  neighbors = 4,
  wrap = TRUE
)

# Display it as a graph
plot(kubo)

# A fun plant model
mod <- camodel(
  transition("plant", "empty", ~ death * ( 1 - (2*q["plant"]-1)^2) ),
  transition("empty", "plant", ~ q["plant"]^2 ),
  all_states = c("empty", "plant"),
  wrap = TRUE,
  neighbors = 4,
  parms = list(death = 0.2496)
)


# Conway's Game of Life
mod <- camodel(
  transition("LIVE", "DEAD", ~ q["LIVE"] < (2/8) | q["LIVE"] > (3/8)),
  transition("DEAD", "LIVE", ~ q["LIVE"] == (3/8)),
  wrap = TRUE,
  neighbors = 8,
  all_states = c("DEAD", "LIVE")
)

# A spiral-generating rock-paper-scissor model
mod <- camodel(
  transition(from = "r", to = "p", ~ q["p"] > 0.25 ),
  transition(from = "p", to = "c", ~ q["c"] > 0.25 ),
  transition(from = "c", to = "r", ~ q["r"] > 0.25 ),
  parms = list(prob = 1),
  wrap = TRUE,
  neighbors = 8
)

# Display the model as a graph
plot(mod)

# Running the above model (see also the help files for the relevant functions)
init <- generate_initmat(mod, c(r = 1/3, p = 1/3, c = 1/3), nrow = 128)
out <- run_camodel(mod, init, times = seq(0, 128))
plot(out)

# Update a model definition using update()
mod <- camodel(
  transition("plant", "empty", ~ m),
  transition("empty", "plant", ~ r * q["plant"] * ( 1 - q["plant"] ) ),
  all_states = c("empty", "plant"),
  wrap = TRUE,
  neighbors = 4,
  parms = list(m = 0.35, r = 0.4)
)

mod_updated <- update(mod, parms = list(m = 0.05, r = 1))
init <- generate_initmat(mod_updated, c(plant = 0.8, empty = 0.2), nrow = 128)
out <- run_camodel(mod_updated, init, times = seq(0, 128))
plot(out)
image(out)

# You can also specify only part of the parameters, the others will be
# kept to their original values
mod_updated <- update(mod, parms = list(m = 0.035))