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
= birth probability.
= spread risk.
= 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:
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,
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,
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,
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.
# 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, 100def 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 inrange(1,nx-1):
for iy inrange(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
breakelse:
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:
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.
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,
At each chronon, a fish moves randomly to one of the adjacent unoccupied squares. If there are no free squares, no movement takes place.
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,
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.
At each chronon, each shark is deprived of a unit of energy.
Upon reaching zero energy, a shark dies.
If a shark moves to a square occupied by a fish, it eats the fish and earns a certain amount of energy.
Once a shark has survived a certain number of chronons it may reproduce in exactly the same way as the fish.
### 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=idself.x, self.y = x, y
self.energy = init_energy
self.fertility_threshold = fertility_threshold
self.fertility =0self.dead =Falseclass 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 inrange(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 inrange(ncreatures):
whileTrue:
x, y =divmod(random.randrange(self.ncells), self.height)
ifnotself.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].idifself.grid[y][x] else0for x inrange(self.width)] for y inrange(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 =Falseif 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 +=2self.grid[yp][xp].dead =Trueself.grid[yp][xp] = EMPTY
moved =TrueexceptIndexError:
# No fish to eat: just move to a vacant cell if possible.passifnot moved:
# Try to move to a vacant celltry:
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 =TrueexceptIndexError:
# Surrounding cells are all full: no movement.
xp, yp = creature.x, creature.y
if creature.energy <0:
# Creature dies.
creature.dead =Trueself.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 =0self.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 inrange(ncreatures):
creature =self.creatures[i]
if creature.dead:
# This creature has been eaten so skip it.continueself.evolve_creature(creature)
# Remove the dead creaturesself.creatures = [creature for creature inself.creatures
ifnot 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 inrange(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')
Extinction of sharks: if initial number of fish is low, sharks starve to death and become extinct. Then fish flourish.
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.
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
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.
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).
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)
steady-states
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)
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)
Exercises
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.
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.
Specify the rules of the cellular automata in enough detail that we can program it.
Describe the emergent phenomena exhibited by your automata. (waves, particles, solitons, oscillators, interfaces, spirals, freezing, clustering, bursting, ...)