# Lattices for space (2 days)

( previous, home, next )

## Cellular automata

The percolation model is a definitely a spatial model but is rather dull since it doesn't change over time. Lots of problems involve some description of how a spatial pattern changes over time. Perhaps the simplest examples of such games are cellular automata.

Cellular automata created by Stanislaw Ulam and John von Neumann in the 1940's at Los Alamos National Lab as simple models of computation Stan and John were working with some of the first large-scale programmable computers, and there was allot of creativity going on trying to figure out how to use them. A cellular automata consists of a space that's been broken up into discrete "cells", and a rule that specifies how each cell is updated from one time step to the next. In their most common forms, the update rule is deterministic (does not include randomness) and the cells are squares arranged in a line or grid that are all updated simultaneously.

## Forest fire model

Browser no support HTML5 Canvas.

= birth probability.

= probability of lighting.

Re-imagine our percolation problem as a fire burning through a forest. The fire is percolating through.

Let's re-frame a bit. Each cell is either empty, contains a patch of trees, or is on fire. The first rule for this cellular automata model:

1. For a tree-patch cell, if neighbor cell is on fire at time $$t$$, the cell is on fire at time $$t+1$$.

Of course, the fire doesn't burn forever! The second rule is,

1. If a cell is on fire at time $$t$$, it is empty at time $$t+1.$$
Assume the probability that a cell is occupied by a tree patch is 0.5:

This is reasonable, short term. We can use data to inform the initial configuration (density of trees) and ask, what's the probability of the fire getting from one side of the woods to the other?

But more long term, we have additional dynamics. For one, trees grow back. But not immediately. Third rule,

1. An empty space fills with a tree patch with probability $$p$$.

Think of tree regrowth as a Poisson process, so $$p$$ is the probability of regrowth per unit time, say per year.

Even if your initial conditions contain fire patches, these three rules will ultimately result in a dense forest (all cells of your CA model will contain trees). But that's not realistic either. There are many things that will control tree density - resource competition, for one. But an important factor is more forest fires! Lightning strikes cause forest fires. Fourth rule,

1. A tree ignites with probability $$f$$ even if no neighbors are on fire.

Again, think of this as a Poisson process, so $$f$$ is the probability of lightning strike in some unit time, say per year.

Take $$p=0.05$$ and $$f=0.001$$:

That's it! That's our forest fire model. An important control parameter is the ratio $$p/f$$, average number of trees that regrow between lightning strikes.

What happens when $$p/f$$ is small? Take $$p=1/365/5$$ and $$f=1/365$$:

What's more reasonable? Simulate and justify to yourself.

[Show code]
 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  # Modified from http://scipython.com/blog/the-forest-fire-model/ import numpy as np import matplotlib.pyplot as plt from matplotlib import animation from matplotlib import colors # Displacements from a cell to its eight nearest neighbours neighbourhood = ((-1,-1), (-1,0), (-1,1), (0,-1), (0, 1), (1,-1), (1,0), (1,1)) # Moore #neighbourhood = ((-1,0), (0,-1), (0, 1),(1,0)) # von Neumann EMPTY, TREE, FIRE = 0, 1, 2 # Colours for visualization: brown for EMPTY, dark green for TREE and orange # for FIRE. Note that for the colormap to work, this list and the bounds list # must be one larger than the number of different values in the array. colors_list = [(0.2,0,0), (0,0.5,0), (1,0,0), 'orange'] cmap = colors.ListedColormap(colors_list) bounds = [0,1,2,3] norm = colors.BoundaryNorm(bounds, cmap.N) # Forest size (number of cells in x and y directions). nx, ny = 100, 100 def iterate(X): """Iterate the forest according to the forest-fire rules.""" # The boundary of the forest is always empty, so only consider cells # indexed from 1 to nx-2, 1 to ny-2 X1 = np.zeros((ny, nx)) for ix in range(1,nx-1): for iy in range(1,ny-1): if X[iy,ix] == EMPTY and np.random.random() <= p: X1[iy,ix] = TREE if X[iy,ix] == TREE: X1[iy,ix] = TREE for dx,dy in neighbourhood: if np.random.random() <= g: if X[iy+dy,ix+dx] == FIRE: X1[iy,ix] = FIRE break else: if np.random.random() <= f: X1[iy,ix] = FIRE return X1 # The animation function: called to produce a frame for each generation. def animate(i): im.set_data(animate.X) animate.X = iterate(animate.X) # The initial fraction of the forest occupied by trees. forest_fraction = 0.45 # Probability of new tree growth per empty cell, and of lightning strike, and of fire spread p, f, g= 1./365./5., 1./365., 1 # Initialize the forest grid. X = np.zeros((ny, nx)) X[1:ny-1, 1:nx-1] = np.random.randint(0, 2, size=(ny-2, nx-2)) X[1:ny-1, 1:nx-1] = np.random.random(size=(ny-2, nx-2)) < forest_fraction X[1,int(nx/2.)] = 2 #fig = plt.figure(figsize=(25/3, 6.25)) #ax = fig.add_subplot(111) #ax.set_axis_off() #im = ax.imshow(X, cmap=cmap, norm=norm)#, interpolation='nearest') #plt.show() # fig = plt.figure(figsize=(25/3, 6.25)) ax = fig.add_subplot(111) ax.set_axis_off() im = ax.imshow(X, cmap=cmap, norm=norm)#, interpolation='nearest') # The animation function: called to produce a frame for each generation. def animate(i): im.set_data(animate.X) animate.X = iterate(animate.X) ## Bind our grid to the identifier X in the animate function's namespace. animate.X = X # Interval between frames (ms). interval = 100 anim = animation.FuncAnimation(fig, animate, interval=interval) plt.show()

Notes:

1. This model can exhibit a fractal frequency-size distribution of clusters, albeit in a physically irrelevant regime, $$f\ll p\ll T$$ where $$T$$ is the burn time of the largest cluster. Need $$f\ll p$$ so that large clusters can form, and $$p\ll T$$ so that no trees pop up alongside a cluster while it burns. But $$T$$ is not trivial to obtain.

2. Replace "tree" with "susceptible", "on fire" with "infected", and "empty" with "recovered", and change interpretations of $$p$$, $$f$$. And there you have a simple CA SIR model!

## Wa-Tor

Sharks and fish wage an ecological war on the toroidal planet Wa-Tor

"Somewhere, in a direction that can only be called recreational at a distance limited only by one's programming prowess, the planet Wa-Tor swims among the stars. It is shaped like a torus, or doughnut, and is entirely covered with water. The two dominant denizens of Wa-Tor are sharks and fish, so called because these are the terrestrial creatures they most closely resemble. The sharks of Wa-Tor eat the fish and the fish of Wa-Tor seem always to be in plentiful supply."

Original article by A.K. Dewdney, in Scientific American (1984) link.

### Toroidal planet

From original article:

Therefore we can simulate the planet Wa-Tor as a flat sheet with periodic boundary conditions on top/bottom, right/left.

### Fish vs sharks

Simulate fish vs shark dynamics on a grid:

All we need now is rules...

### Rules

"Time passes in discrete jumps, which I shall call chronons. During each chronon a fish or shark may move north, east, south or west to an adjacent point, provided the point is not already occupied by a member of its own species. A random-number generator makes the actual choice. For a fish the choice is simple: select one unoccupied adjacent point at random and move there. If all four adjacent points are occupied, the fish does not move. Since hunting for fish takes priority over mere movement, the rules for a shark are more complicated: from the adjacent points occupied by fish, select one at random, move there and devour the fish. If no fish are in the neighborhood, the shark moves just as a fish does, avoiding its fellow sharks."

#### For the fish,

1. At each chronon, a fish moves randomly to one of the adjacent unoccupied squares. If there are no free squares, no movement takes place.

2. Once a fish has survived a certain number of chronons it may reproduce. This is done as it moves to a neighbouring square, leaving behind a new fish in its old position. Its reproduction time is also reset to zero.

#### For the sharks,

1. At each chronon, a shark moves randomly to an adjacent square occupied by a fish. If there are none, the shark moves to a random adjacent unoccupied square. If there are no free squares, no movement takes place.

2. At each chronon, each shark is deprived of a unit of energy.

3. Upon reaching zero energy, a shark dies.

4. If a shark moves to a square occupied by a fish, it eats the fish and earns a certain amount of energy.

5. Once a shark has survived a certain number of chronons it may reproduce in exactly the same way as the fish.

### Simulations

[Show code]
 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 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260  ### Modified from http://scipython.com/blog/wa-tor-world/ import random import matplotlib.pyplot as plt from matplotlib.colors import LinearSegmentedColormap from matplotlib import animation from matplotlib import colors EMPTY = 0 FISH = 1 SHARK = 2 # Colour the cells for the above states in this order: colors = ['lightblue', 'darkorange', 'grey'] n_bin = 3 cm = LinearSegmentedColormap.from_list( 'wator_cmap', colors, N=n_bin) # Run the simulation for MAX_CHRONONS chronons (time intervals). MAX_CHRONONS = 400 # Save every SAVE_EVERYth chronon iteration. SAVE_EVERY = 5 # PRNG seed. SEED = 20 random.seed(SEED) initial_energies = {FISH: 20, SHARK: 3} fertility_thresholds = {FISH: 4, SHARK: 12} class Creature(): """A sea creature living in Wa-Tor world.""" def __init__(self, id, x, y, init_energy, fertility_threshold): """Initialize the creature. id is an integer identifying the creature. x, y is the creature's position in the Wa-Tor world grid. init_energy is the creature's initial energy: this decreases by 1 each time the creature moves and if it reaches 0 the creature dies. fertility_threshold: each chronon, the creature's fertility increases by 1. When it reaches fertility_threshold, the creature reproduces. """ self.id = id self.x, self.y = x, y self.energy = init_energy self.fertility_threshold = fertility_threshold self.fertility = 0 self.dead = False class World(): """The Wa-Tor world.""" def __init__(self, width=100, height=100): """Initialize (but don't populate) the Wa-Tor world.""" self.width, self.height = width, height self.ncells = width * height self.grid = [[EMPTY]*width for y in range(height)] self.creatures = [] def spawn_creature(self, creature_id, x, y): """Spawn a creature of type ID creature_id at location x,y.""" creature = Creature(creature_id, x, y, initial_energies[creature_id], fertility_thresholds[creature_id]) self.creatures.append(creature) self.grid[y][x] = creature def populate_world(self, nfish=320, nsharks=107): """Populate the Wa-Tor world with fish and sharks.""" self.nfish, self.nsharks = nfish, nsharks def place_creatures(ncreatures, creature_id): """Place ncreatures of type ID creature_id in the Wa-Tor world.""" for i in range(ncreatures): while True: x, y = divmod(random.randrange(self.ncells), self.height) if not self.grid[y][x]: self.spawn_creature(creature_id, x, y) break place_creatures(self.nfish, FISH) place_creatures(self.nsharks, SHARK) def get_world_image_array(self): """Return a 2D array of creature type IDs from the world grid.""" return [[self.grid[y][x].id if self.grid[y][x] else 0 for x in range(self.width)] for y in range(self.height)] def get_world_image(self): """Create a Matplotlib figure plotting the world.""" im = self.get_world_image_array() fig = plt.figure(figsize=(8.3333, 6.25), dpi=72) ax = fig.add_subplot(111) ax.imshow(im, interpolation='nearest', cmap=cm) # Remove ticks, border, axis frame, etc ax.set_xticks([]) ax.set_yticks([]) ax.axis('off') return fig def show_world(self): """Show the world as a Matplotlib image.""" fig = self.get_world_image() plt.show() # plt.close(fig) def save_world(self, filename): """Save a Matplotlib image of the world as filename.""" fig = self.get_world_image() # NB Ensure there's no padding around the image plot plt.savefig(filename, dpi=72, bbox_inches='tight', pad_inches=0) plt.close(fig) def get_neighbours(self, x, y): """Return a dictionary of the contents of cells neighbouring (x,y). The dictionary is keyed by the neighbour cell's position and contains either EMPTY or the instance of the creature occupying that cell. """ neighbours = {} for dx, dy in ((0,-1), (1,0), (0,1), (-1,0)): xp, yp = (x+dx) % self.width, (y+dy) % self.height neighbours[xp,yp] = self.grid[yp][xp] return neighbours def evolve_creature(self, creature): """Evolve a given creature forward in time by one chronon.""" neighbours = self.get_neighbours(creature.x, creature.y) creature.fertility += 1 moved = False if creature.id == SHARK: try: # Try to pick a random fish to eat. xp, yp = random.choice([pos for pos in neighbours if neighbours[pos]!=EMPTY and neighbours[pos].id==FISH]) # Eat the fish. Yum yum. creature.energy += 2 self.grid[yp][xp].dead = True self.grid[yp][xp] = EMPTY moved = True except IndexError: # No fish to eat: just move to a vacant cell if possible. pass if not moved: # Try to move to a vacant cell try: xp, yp = random.choice([pos for pos in neighbours if neighbours[pos]==EMPTY]) if creature.id != FISH: # The shark's energy decreases by one unit when it moves. creature.energy -= 1 moved = True except IndexError: # Surrounding cells are all full: no movement. xp, yp = creature.x, creature.y if creature.energy < 0: # Creature dies. creature.dead = True self.grid[creature.y][creature.x] = EMPTY elif moved: # Remember the creature's old position. x, y = creature.x, creature.y # Set new position creature.x, creature.y = xp, yp self.grid[yp][xp] = creature if creature.fertility >= creature.fertility_threshold: # Spawn a new creature and reset fertility. creature.fertility = 0 self.spawn_creature(creature.id, x, y) else: # Leave the old cell vacant. self.grid[y][x] = EMPTY def evolve_world(self): """Evolve the Wa-Tor world forward in time by one chronon.""" # Shuffle the creatures grid so that we don't always evolve the same # creatures first. random.shuffle(self.creatures) # NB The self.creatures list is going to grow as new creatures are # spawned, so loop over indices into the list as it stands now. ncreatures = len(self.creatures) for i in range(ncreatures): creature = self.creatures[i] if creature.dead: # This creature has been eaten so skip it. continue self.evolve_creature(creature) # Remove the dead creatures self.creatures = [creature for creature in self.creatures if not creature.dead] # advance the world 1 time step tots = [] def advance(): world.evolve_world() X = world.get_world_image_array() tots.append([sum([ i.count(j) for i in X]) for j in range(3)]) return X # The animation function: called to produce a frame for each generation. def animate(i): if i<=1000: im.set_data(animate.X) animate.X = advance() # Initialize world = World() world.populate_world() ## Animation time. First plot IC X = world.get_world_image_array() fig = plt.figure(figsize=(25/3, 6.25)) ax = fig.add_subplot(111) ax.set_axis_off() im = ax.imshow(X, cmap=cm) #, interpolation='nearest' ## Bind our grid to the identifier X in the animate function's namespace. animate.X = X # Interval between frames (ms). interval = 100 anim = animation.FuncAnimation(fig, animate, interval=interval) plt.show() #anim.save('WaTor.mp4',fps=15) # only works outside of canopy ######################################################## # now plot the predator-prey dynamics tots=array(tots) plot(tots[:,1],'darkorange',linewidth=3,label='Fish') plot(tots[:,2],'grey',linewidth=3,label='Sharks') xlabel('Time (chrons)',fontsize=24) ylabel('Pop. sizes',fontsize=24) xticks(fontsize=20) yticks(fontsize=20) axis([0,1000,0,10000]) legend(loc=1,fontsize=20) plt.tight_layout() savefig('WaTor_TimeSeries.png')

#### Possible outcomes

Light blue = empty water, grey = sharks, orange = fish.

1. Extinction of sharks: if initial number of fish is low, sharks starve to death and become extinct. Then fish flourish.
2. Extinction of both species: if the number of sharks is so great that they consume all the fish, i.e., drive fish to extinction, then starve, going extinct themselves.
3. Coexistence of sharks and fish.

Do these outcomes sound familiar?...

#### Focus on coexistence case

Sharks and fish populations on Wa-Tor, plotted over time:

Note the oscillations... remind you of anything?

We have simulated a spatial predator-prey model. The spatial modeling has revealed complexity inaccessible in our previous, ODE model. We see different patterns forming, the hints of traveling waves, spiral waves, etc. Go back and watch the movie, run your own simulations, and see for yourself. Thus we gain a deeper understanding of predator-prey systems.

Bifurcations, or some CA-analogy, which would help us understand basins of attractions for these outcomes, are a recent and ongoing area of research.

## Heart Rhythms

Lattice--Boltzmann models of electrical signal transmission in the heart

Rules

• A cell that is currently firing goes dark (refractory) in the next step and stays that way for 10 steps.

• Then a cell becomes blue (fresh) and stays that way unless it is triggered by an active firing neighbor.

• If a neighbor of a fresh cell is red (actively firing), then there is an 85 % chance that fresh cell will actively fire in the next step.

Your browser does not support HTML5 Canvas.

= Probability of firing.

## Artificial life

One of the best known cellular automata models is Conway's "game of life". The game of life is played on a square lattice using a Moore neighborhood. Each site takes one of two states -- occupied or empty. The transition rules for the game of life are totalistic -- they depend only on the number of neighbors, not the specific configuration of neighbors.

The rules can be loosely thought of as modelling the life of single-celled organisms in a petri dish.

• Use a Moore neighborhood on a rectangular lattice with wrapped boundaries to make a torus.
• Any occupied site with 2 or 3 occupied neighbours remains occupied.
• Any occupied site with fewer than 2 occupied neighbors becomes empty (under-population).
• Any occupied site with more than 3 occupied neighbours empties, (overcrowding).
• Any empty site with exactly three occupied neighbours becomes occupied (reproduction).

These rules create a menagerie of phenomena.

### 1-d automata

Today, cellular automata are most often associated with the work of Stephen Wolfram, the physicist who created Mathematica. Wolfram conducted a systematic study and classification of the simplest interesting cellular automata -- binary line automata with only nearest neighbor interactions.

Rule enumeration method

compact state space -> must eventually repeat. (what will the period be?)

Example: ??

Classifications

• period (all will have finite period)
• attractors independent of initial condition?
• attractors fix parts of initial condition (error-correcting)
• amphichiral oscillators (period 2, 4, ...)
• shifting patterns (period corresponds to domain harmonics)
• long-period acyclics (automata chaos)

## Schelling's suburbs

Visual implementation at Case2014

Model of how technology and network effects change the world in irreversible ways. Two competing systems, each with positive minimum capital thresholds, people start switching for one reason, others are forced to switch, get locked in, old established system collapses, and nobody can switch back any more.

## Voter model

• 8-cell Moore neighborhoods on a square lattice.
• Every cell represents a voter with one of two opinions on a question.
• At each step, one voter is picked at random, and that voter changes their opinion to match that of one of their neighbors.

= # steps between refreshes (requires stop/start)

Browser no support HTML5 Canvas.

# Exercises

1. Chirality. One of the interesting things about Conway's game of life and related totalistic cellular automata is that the rules are not just translationally and rotationally symmetric, but also achiral -- lattice animals that curve to the right get transformed the same way as if they were flipped so they curved to the left. But what if the rules were symmetric but chiral, so there is a natural curvature leads to different shape evolutions -- would the dynamics have any new weird properties? Chirality turns out to be a important thing -- in through the looking glass, it's commented that milk tastes different on the other side -- because of chirality. Make a cellular automata rule that is rotationally symmetric but achiral, simulate the rule, and discuss.

2. Search the web for a new cellular-automata model that we have not already discussed in class. Your rule should be different from your classmates.

1. Specify the rules of the cellular automata in enough detail that we can program it.

2. Describe the emergent phenomena exhibited by your automata. (waves, particles, solitons, oscillators, interfaces, spirals, freezing, clustering, bursting, ...)

( previous, home, next )