# Conway's Game of Life

## Set-up

The first step of running a Julia program is to load the external libraries one will be using. We do this with a `using`

statement.

`using AlgebraicABMs, Catlab, AlgebraicRewriting, Random, Test`

## Schema

Defining an schema is stating what data is required to specify a state of the simulation at some point in time. In AlgebraicJulia, this is done via declaring a `Presentation`

, i.e. a database schema. Objects (`Ob`

, or tables) are types of entities. Homs (`Hom`

, or foreign keys) are functional relationships between the aforementioned entities. AttrTypes are placeholders for Julia types, which are assigned to objects via attributes (`Attr`

).

The schema below *extends* the schema for symmetric graphs, which consists in two tables (`E`

and `V`

, for edges and vertices), a hom `inv`

which pairs each edge with its symmetric dual edge, and homs (`src`

, `tgt`

) which relate the edges to the vertices. We consider each vertex to be a cell in the Game of Life. This is a generalization of the game, since it is not enforced that every cell has eight neighbors. However the example we run at the end will be a regular grid of vertices.

Everything in the `@present`

presentation below simply adds to the schema of symmetric graphs, indicated by the subtype operator `<: SchSymmetricGraph`

. We need one more piece of information to specify a state of the world in the game of life: which cells are alive? There are a few ways one could distinguish the living cells from the dead ones. Here, we add a new table, `Life`

which can be though of as a set of "life tokens". A function, `live`

, assigns each such token to a cell, designating it to be alive. Thus the subset of living cells is the image of the `live`

function.

To summarize, a state of the world in the Game of Life is a set of cells (i.e. vertices), a set of edges (with `src`

/`tgt`

functions) which encode which cells are near each other, and a distinguished subset of cells given by the `live`

function which marks which cells are alive. This information is summarized in the following graphical depiction of the schema.

```
@present SchLifeGraph <: SchSymmetricGraph begin
Life::Ob
live::Hom(Life,V)
end
to_graphviz(SchLifeGraph) # visualize the schema
```

If `SchLifeGraph`

is the piece of data which describes what it means to be a world state, we need a Julia datatype whose values are world staets.

The `@acset_type`

macro takes the specification of a datatype (i.e. a schema), and it generates an efficient Julia datatype, which is essentially the type of in-memory databases which have their schema given by `SchLifeGraph`

.

This line defines `LifeState`

to be this new Julia type, and futhermore it states that it satisfies the `AbstractSymmetricGraph`

interface: Catlab already has a lot of generic code to work with SymmetricGraph-like things, so now we are free to use it with LifeState.

`@acset_type LifeState(SchLifeGraph) <: AbstractSymmetricGraph;`

Before we start defining the dynamics of our Game of Life, it'll be helpful to visualize states of the world. One thing that `SchLifeGraph`

did not have was (x,y) coordinate information. Because we were representing proximity of vertices by the prescence or absence of an edge connecting them, these coordinates are irrelevant to the model dynamics.

However, the coordinates are useful for visualization, so the below code creates a new schema which extends the previous one to add coordinates.

```
@present SchLifeCoords <: SchLifeGraph begin
Coords::AttrType
coords::Attr(V, Coords)
end
to_graphviz(SchLifeCoords)
```

Likewise for this new schema, we need to use `@acset_type`

to turn the description of the required data of a world state (the schema) into an actual Julia datatype. One difference here is that there is now an attribute present, which means we need to indicate what Julia type it should be.

Our new Julia type is `LifeStateCoords`

, and we pick `Tuple{Int,Int}`

as our representation of (x,y) coordinates.

`@acset_type LifeStateCoords(SchLifeCoords){Tuple{Int,Int}} <: AbstractSymmetricGraph;`

The following (hidden) code visualizes LifeStates and uses coordinates if the input is a `LifeStateCoords`

, but otherwise it places the vertices in an arbitrary location.

`function view_life(X::Union{LifeState, LifeStateCoords}, pth=tempname())`

This (hidden) helper function creates an ordinary, square grid and initializes cells as alive or dead via a boolean-valued input matrix. Like above, one does not need to understand this in order to understand the Game of Life model - it's just a convenient way to generate some initial states to run the model on.

`function make_grid(curr::AbstractMatrix)`

Another helper function, which creates a game state on a square n × n grid. If the keyword `random`

is `true`

(which it is, by default), then cells have a 50% chance of being alive.

Note: red color means that a cell is not in the image of the `live`

function.

```
make_grid(n::Int, random=true) =
make_grid((random ? rand : zeros)(Bool, (n, n)));
view_life(make_grid(3)) # For example, visualize a random 3x3 grid
```

We've now constructed two schemas and two datatypes. No doubt they are related to each other, and AlgebraicJulia gives us the ability to take advantage of that relationship in order to automatically values of `LifeStateCoords`

into values of `LifeState`

(it does this is the obvious way: it throws away the coordinate information). More interestingly, it gives us a way of converting a `LifeState`

into a `LifeStateCoords`

: we give "generic" coordinates to each of the vertices.

If we call `AddCoords`

on something built out of `LifeState`

instances, it will make the corresponding thing built out of `LifeStateCoords`

instances.

`AddCoords = Migrate(SchLifeGraph, LifeState, SchLifeCoords, LifeStateCoords; delta=false);`

## Initial state for the model

With our `make_grid`

function, it's easy to create an initial state.

```
G = make_grid([0 1 0;
1 1 1;
0 1 0])
view_life(G)
```

We'll use this world-state, called `G`

throughout the model-building process to confirm our model building blocks work the way they're supposed to work.

## Building blocks for the Game of Life model

We are now ready to build our model, which will have three rules: underpopulation, overpopulation, and birth. When we write these rules down, it will be very succinct due to the hard work of this section: laying out the building blocks.

Our first building block is an instance of `LifeState`

. Now, we don't want to think of this as a particular state of the game of life, but rather we want to think of it as a *pattern* for which we could look for *matches* inside a real game state.

This pattern is very simple: it consists of a single vertex. A *match* for this pattern in some world state `X`

is just a choice of a vertex in `X`

. The matched vertex need not be dead just because our pattern vertex is not marked as alive. So the meaning of `Cell`

is "a dead or alive cell", despite the fact that, when we visualize it, it shows up as red (because it hasn't been marked explicitly to be alive).

```
const Cell = LifeState(1) # this means: a graph with one vertex, nothing else
view_life(Cell) # red color means this cell is not in image of `life` function.
```

Let's confirm that this works. Our fancy word for "pattern match" is "homomorphism" ("morphism", for short). To look for morphisms from `Cell`

into `G`

is to ask "find me all the cells!". The Catlab function `homomorphisms`

exhaustively finds all such answers to a query like this. When we run it below, we ought see there are 9 answers to this query, as there are 9 vertices in the 3 × 3 grid.

`length(homomorphisms(AddCoords(Cell), G))`

`9`

Our next building block is also a `LifeState`

instance that we think of as a pattern. It is a cell but, furthermore, it is a cell which has been marked as alive. This means it consists in one vertex, one 'life token', and the function `live`

which maps `Life₁`

to `V₁`

.

```
const LiveCell = @acset LifeState begin V=1; Life=1; live=1 end
view_life(LiveCell) # green color means this cell is in image of `life` function
```

Let's confirm that this works. If we look for morphisms from `LiveCell`

into `G`

, meaning "find me all the live cells!", we will see there are 5 answers to this query

`length(homomorphisms(AddCoords(LiveCell), G))`

`5`

Now we will see our first morphism that is not between some pattern and some world state but, rather, is *between* our the two previous patterns. A morphism from `Cell`

is an answer to the query "give me a Cell!" and there is exactly one possible answer to that query in `LiveCell`

, thus we can use `homomorphism`

, which returns an arbitrary morphism from the set of morphisms computed by `homomorphisms`

.

`to_life`

is an important morphism to understand as it is the fundamental "engine" of the Game of Life. When read 'forwards', it expresses the idea of taking a vertex and adding a Life token to it. When read 'backwards', it expresses the idea of removing a Life token from a vertex. This is all we actually do when simulating the Game of Life: we toggle cells between being marked as alive vs no longer being marked as alive.

`const to_life = homomorphism(Cell, LiveCell);`

The remaining challenge is to control the circumstances under which we perform this fundamental action of toggling life status. This is done via a technique called "Positive (resp. negative) application conditions". An application condition is a way of embedding a small pattern in a larger pattern (a context) with the assertion that, for a specific match of the small pattern in a world state, that this context is required (resp. forbidden) to hold.

Here the small pattern, the large pattern, and the world state will all be `LifeState`

instances. And the the 'embedding' of the small pattern into the large pattern as well as the match of the small pattern in the world will both be morphisms between `LifeState`

instances.

The following lines of code provide shorthand names for positive and negative applications. The use of the keyword `monic`

can be ignored if this is your first time looking at AlgebraicJulia code. If not, the keyword indicates that the morphism (from the larger context to the world state, used to determine whether the application condition is satisfied) is required to be monic.

```
PAC(m) = AppCond(m; monic=true) # Positive Application condition
NAC(m) = AppCond(m, false; monic=true); # Negative Application condition
```

The following line provides a shorthand for constructing the fundamental building block of an ABM, which is an `ABMRule`

. Such a rule requires a name, an instance of a rewrite `Rule`

object from `AlgebraicRewriting`

, as well as a timer which indicates the frequency at which the rule fires. All three of our rules will fire at every tick in time.

```
TickRule(name, I_L, I_R; ac) = # Rule which fires on 1.0, 2.0, ...
ABMRule(name, Rule(I_L, I_R; ac), DiscreteHazard(1));
```

The context for determining whether the rules of underpopulation, overpopulation, or birth apply are always framed in terms of "how many living neighbors do you have?". Thus, it is helpful to have function which creates a context of `n`

living neighbors which we can use with various values of `n`

.

Providing a context for a pattern is more than just giving an instance of `LifeState`

; we also have to say how the pattern relates to the context. So our pattern here is either a vertex (if `alive=false`

) or a vertex marked as alive (`alive=true`

, which is the default value). This morphism data is what allows us to pick out the *center* vertex as the relevant starting point when checking if a particular match (which is a choice of a vertex in the game state) actually exists with that context of `n`

living living neighbors.

When we visualize the result of this function below, we're only visualizing the *codomain* of the morphism (i.e. what the morphism is pointing to).

```
function living_neighbors(n::Int; alive=true)::ACSetTransformation
view_life(codom(living_neighbors(3; alive=false)))
```

## Create model by defining update rules

We now have our building blocks in place and can build the three rules.

A rewrite rule is given by a pair of morphisms, `L ← I → R`

. The first morphism is 'read backwards' (thought of as deletion), and the second one is 'read forwards' (thought of as addition). When one's rule does pure addition, the first morphism is an identity map. Likewise when one's rule does pure deletion, the second morphism is an identity map. So rules `underpop`

and `overpop`

will have their rule given by first `to_life`

and then `id(Cell)`

, whereas the rule `birth`

will first have `id(Cell)`

and then `to_life`

. The real interesting part of these constructions then is their application conditions (given via a keyword argument `ac`

).

Let's start with underpopulation: a cell dies due to underpopulation if it has fewer than two living neighbors. So the application condition is a *negative* application: if we notice a cell has the context of two living neighbors, we should **not** apply the `underpop`

rule.

```
underpop =
TickRule(:Underpop, to_life, id(Cell); ac=[NAC(living_neighbors(2))]);
```

A cell dies due to overpopulation if it has over 3 living neighbors. So here our condition of applying the death rewrite (i.e. reading the morphism `to_life`

backwards) is a *positive* application condition. We cannot kill a cell unless, furthermore, that cell is in a context where it has (at least) four living neighbors.

```
overpop =
TickRule(:Overpop, to_life, id(Cell); ac=[PAC(living_neighbors(4))]);
```

The rule for birth is that one must have a dead cell with exactly three living neighbors. Our pattern here is simple a vertex, which could match any cell in the game state! So we need to use multiple application conditions to make sure the cell has at least three neighbors, fewer than four neighbors, and doesn't already have a life token. These three constraints are respectively the three elements of the list below.

```
birth = TickRule(:Birth, id(Cell), to_life;
ac=[PAC(living_neighbors(3; alive=false)),
NAC(living_neighbors(4; alive=false)),
NAC(to_life)]); # this rule does NOT apply if cell is alive
```

We can now create the model: an ABM is constituted by its transition rules

`GoL = ABM([underpop, overpop, birth]);`

We wrote an ABM for `LifeState`

, but we want to apply it to `G`

as its initial state. `G`

is not a `LifeState`

, but rather a `LifeStateCoords`

! Luckily, we defined `AddCoords`

earlier which knows how to migrate things in the language of `LifeState`

into the language of `LifeStateCoords`

.

`GoL_coords = AddCoords(GoL);`

## Checking our rules make sense

Let's remind ourselves what G looks like

`view_life(G)`

This code below counts how many matches a rule has, application conditions considered. It returns the coordinates of those matches.

```
match_coords(f::ACSetTransformation) = codom(f)[f[:V](1), :coords]
match_coords(rule::ABMRule, X) = match_coords.(get_matches(AddCoords(rule), X));
```

Let's calculate which cells will die from underpopulation in the first time step:

`match_coords(underpop, G)`

`Any[]`

This is right, there are no living cells which have fewer than two living neighbors.

Now, we see that the center cell dies from overpopulation:

`match_coords(overpop, G)`

```
1-element Vector{Tuple{Int64, Int64}}:
(2, 2)
```

Below are the coordinates of the four cells that will come to life:

`match_coords(birth, G)`

```
4-element Vector{Tuple{Int64, Int64}}:
(1, 1)
(3, 1)
(1, 3)
(3, 3)
```

## Running the ABM

This is easy! Pass in our model and our initial state. We optionally could limit the run via some maximum number of steps or time, but this one will achieve steady state within 3 time steps.

```
res = run!(GoL_coords, G);
@test length(res) == 13
```

`Test Passed`

## View results

We use `view_life`

to generate a bunch of images capturing the trajectory.

`imgs = view(res, view_life);`

The first image is our starting point.

`imgs[1]`

The next steps (all taking place at `t=1.0`

) are the four births and one overpopulation death.

`imgs[2]`

`imgs[3]`

`imgs[4]`

`imgs[5]`

`imgs[6]`

At this point, we move to `t=2.0`

. There are now four cells which are overpopulated.

`imgs[7]`

`imgs[8]`

`imgs[9]`

`imgs[10]`

In the last step, the remaining four living cells perish due to underpopulation.

`imgs[11]`

`imgs[12]`

`imgs[13]`

`imgs[14]`

## Bonus: stochastic game of life

Rather than having all cells update in lockstep, we could change the probability of firing from a Dirac delta distribution at t=1 to an exponential distribution, where the expected value is firing at t=1. This means cells will update one at a time, as it is almost impossible for two events to occur at the same time. The change involved for this is simply replacing the `DiscreteHazard`

of `TickRule`

with a `ContinuousHazard`

.

```
continuous_abm = ABM([ABMRule(r.name, r.rule, ContinuousHazard(1)) for r in GoL_coords.rules])
res = run!(continuous_abm, make_grid(4); maxevent=3)
imgs = view(res, view_life);
```

```
┌ Debug: Step 0: Event Underpop | Fired @ t = 0.21 (3 queued)
└ @ AlgebraicABMs.ABMs ~/work/AlgebraicABMs.jl/AlgebraicABMs.jl/src/ABMs.jl:521
┌ Debug: Step 1: Event Underpop | Fired @ t = 0.45 (2 queued)
└ @ AlgebraicABMs.ABMs ~/work/AlgebraicABMs.jl/AlgebraicABMs.jl/src/ABMs.jl:521
┌ Debug: Step 2: Event Underpop | Fired @ t = 1.11 (1 queued)
└ @ AlgebraicABMs.ABMs ~/work/AlgebraicABMs.jl/AlgebraicABMs.jl/src/ABMs.jl:521
```

Here is our starting point.

`imgs[1]`

Let's look at the next few steps.

`imgs[2]`

Then

`imgs[3]`

Then

`imgs[4]`

*This page was generated using Literate.jl.*