-
Notifications
You must be signed in to change notification settings - Fork 5
/
Copy path04-project-1DCA.jmd
158 lines (113 loc) · 6.34 KB
/
04-project-1DCA.jmd
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
---
title : "Project: 1D elementary cellular automata"
author : Michiel Stock, Bram De Jaegher, Daan Van Hauwermeiren
date: 18 December 2019
---
In this exercise, we simulate elementary one-dimensional cellular automata (CA). It is a simple application, which nevertheless illustrates many of the Julia programming concepts as it has all aspects of an individual-based model.
# What are elementary cellular automata?
Cellular automata are simple dynamical models, discrete in time, space, and state. In elementary cellular automata, the states only take the values 0 or 1. The state of the cell at time $t+1$ is solely determined by the state of that cell at time $t$ and the states of its neighbors. Since these three cells can jointly only take eight unique values, i.e.,
```julia
[(l, s, r) for l in 0:1 for s in 0:1 for r in 0:1]
```
This means that if we assign one state to each combination, we have $2^8=256$ unique elementary cellular automata. Naturally, many of them are equivalent up to some symmetry, but there is quite some diversity in behavior among these different classes.
An integer between 0 and 255 solely defines every elementary cellular automaton. The binary representation of this number determines the state updates. This is best explained by showing some examples.
So starting from an initial state vector, we can represent the evolution of the states by a space-time diagram.
In implementing an elementary CA, we use cyclic boundary conditions. Here, we wrap the state vector around a circle: the right neighbor of the last cell is again the first cell (and *vice versa*).
# A basic implementation
## Update rule
We start by implementing the updating rule from a single state to the next. This can be implemented in a variety of ways, but let us provide one that is reasonably concise and efficient. It will also nicely illustrate some aspects of Julia's type system.
Updating a cell depends only on the state of the cell and its two neighbors. Let us use three Boolean variables `(l, s, r)` to denote them. Here `s` is the state of the cell of interest, and `l` and `r` are the states of its left and right neighbor, respectively. An integer between 0 and 255 can represent the rule, hence an 8-bit number. The rule for updating a state is implemented below.
```julia
nextstate(l::Bool, s::Bool, r::Bool, rule::UInt8) = bitstring(rule)[8-(4r+2s+1r)] == '1';
```
Here, we have cleverly made use of the function `bitstring`, which transforms an 8-bit number in the corresponding bitstring.
```julia
bitstring(UInt8(42))
```
The three Boolean variables are transformed into the correct index of this string. We might save ourselves the bother of ensuring that the argument `rule` is of the type `UInt8` using dispatch.
```julia
nextstate(l::Bool, s::Bool, r::Bool, rule::Integer) = nextstate(l, s, r, UInt8(rule));
```
> **Assignment**: Create an 8 by 3 matrix `S` containing all possible states of three adjacent cells. Complete the function `nextstates`, which returns a vector containing the next state of the middle cell.
```julia;eval=false
rule = ...
```
```julia;eval=false
S = ...
```
```julia;eval=false
function nextstates(S, rule)
...
return states
end
```
## Simulating a 1D CA model
We will represent the states of our system using a binary vector (type is `BitVector`) `x`.
Efficient functions often work in-place, because memory is already allocated. We will construct a function `next!`, which takes the current state vector `x` and a vector to store the next states in a binary vector `xnew` for a given rule. This function uses cyclic boundaries.
> **Assignment**: Complete the function `next!`.
```julia;eval=false
"""
next!(xnew, x, rule)
Updates a state vector `x` according to `rule` and stores the result in `xnew`.
"""
function next!(xnew, x, rule)
...
return xnew
end
```
> **Optional assignment**: Make a simular function `next`, which only takes as arguments `x` and `rule`, i.e. this function allocates a new vector for the result.
Now that we can simulate a single time step, let us write a function to simulate for several time steps.
```julia;eval=false
"""
simulate(x₀, rule; nsteps::Integer=100)
Simulate `nsteps` time steps according to `rule` with `x₀` as the initial condition.
Returns a matrix X, where the rows are the state vectors at different time steps.
"""
function simulate(x₀, rule; nsteps::Integer=100)
...
return X
end
```
> **Assignments**
> 1. Complete the function `simulate`.
> 2. Choose a number for the rule. Simulate a elementary CA of size 101 for 100 timesteps. Start from a vector with all elements false, except for the middle cell.
> 3. Visualize the result. This can be done either using `heatmap(1X)` (multiplying by 1 turns the Booleans into numerical values) or the function `printca` given below.
```julia; eval=false
x₀ = ...
```
```julia
function printca(X::Matrix{Bool})
n, m = size(X)
for i in 1:n
println(reduce(*, (s-> s ? "*" : "-").(X[i,:])))
end
end
```
# Type-based dispatch for rules
In the previous section, we developed a working implementation of an elementary cellular automaton. Let us rewrite the code, but using type-based dispatch. We implement an abstract type for general CA rules. Then we have a structure `DetCARule` to contain the above elementary CA rules.
```julia
abstract type CARule end
struct DetCARule <: CARule
rulenr::UInt8
function DetCARule(r::Integer) # custom constructor
new(UInt8(r))
end
end
```
> **Assignment** Complete the code for the type-based dispatch of cellular automata. Test it using the same rule as before (or a different one, we don't care).
```julia;eval=false
nextstate(l::Bool, s::Bool, r::Bool, rule::DetCARule) = ...
function next!(xnew, x, rule::CARule)
...
return xnew
end
function simulate(x₀, rule::CARule; nsteps::Integer=100)
...
return X
end
```
> **Assignment** Make a new concrete type `StochCARule` to model stochastic elementary cellular automata. This type has an additional field `pflip`, which is the probability of switching the state after applying the rule. Make a custom `nextstate` to deal with this and simulate again correctly.
```julia
```