I made a random-walk predator-prey simulation that focuses on individual animals instead of the (maybe) more common array-based approach. I'd like to hear your opinion about this: how could the concept and the code structure be improved, and how good is my implementation?

Rules in the model:

  • There are two species competing on a rectangular grid: rabbits and foxes.
  • Initially an exact number of them is spawned on random distinct locations.
  • In each step, each of them moves 1 step to a neighbour grid (no diagonal movements).
  • If a rabbit and a fox steps on the same spot, the fox eats the rabbit, so the rabbit dies.
  • Each animal starts with a given amount of energy, and each timestep they lose 1 energy. If an animal runs out of energy, it dies.
  • in each step, there's a small chance for each point that a newborn rabbit or fox spawns there, with fresh energy. If a (newborn/not) rabbit and a (newborn/not) fox occupies the same point, the fox eats the rabbit instantly.

And, of course, there's a lot of possibility to play around: change parameters like energy or spawning probability, initial values, rules of interaction (recharge energy when feeding?) or spawning (spawn only if two of the same species on the same place? after feeding)... even adding a new species is fairly simple.

The code:

# -*- coding: utf-8 -*-
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import animation as animation
from scipy import integrate
from progressbar import progressbar as prbar        # (use pip/conda install progressbar2, or rewrite line 116.)
from copy import copy

FOX = 1

UP = 0
DOWN = 1
LEFT = 2
STAY = 4

# initial number of rabbits and foxes
nrabbits = 10
nfoxes = 5

# size of the grid
gridxsize = 30
gridysize = 30

# energy of a freshly spawned rabbit/fox
rabben = 20
foxen = 20

# chance of a new fox/rabbit being spawned at a gridpoint on a step
rabbit_newborn_chance = 0.01
fox_newborn_chance = 0.01

# number of steps to simulate
steps = 200

class Animal(object):
    Tracks the animal's position, energy, species (rabbit/fox) and state (live/dead).

    def __init__(self, x0, y0, init_energy, species):
        self.x = x0
        self.y = y0
        self.energy = init_energy
        self.species = species
        self.isDead = False

    def interact(self, other):
        Interact with another animal:
            - If they're from the same species, ignore each other.
            - Fox eats rabbit.
        if self.species == RABBIT and other.species == FOX:

        elif self.species == FOX and other.species == RABBIT:

    def die(self):
        self.isDead = True

    def move(self, direction):
        """Move a step on the grid. Each step consumes 1 energy; if no energy left, die.
        If hitting the bounds of the grid, "bounde back", step to the opposite direction insetad.

            direction {int} -- direction to move: UP: 0, DOWN: 1, LEFT: 2, RIGHT: 3, STAY: 4
        self.energy -= 1

        if direction == LEFT:
            self.x += 1 if self.x > 0 else -1   #"bounce back"
        if direction == RIGHT:
            self.x -= 1 if self.x < gridxsize-1 else -1
        if direction == UP:
            self.y += 1 if self.y < gridysize-1 else -1
        if direction == DOWN:
            self.y -= 1 if self.y > 0 else -1
        if direction == STAY:

        if self.energy <= 0:
            self.die()          #R.I.P.

animals = []        # this will contain all animals on the grid

# all possible coordinate pair (following https://stackoverflow.com/a/11144716/5099168)
xcoords = np.arange(gridxsize)
ycoords = np.arange(gridysize)
coords = np.transpose([np.tile(xcoords, len(ycoords)), np.repeat(ycoords, len(xcoords))])

# populate grid randomly, unique coordinates for all animals
randcoords = np.random.permutation(coords)
rabbitcoords = randcoords[:nrabbits]
foxcoords = randcoords[nrabbits:(nrabbits + nfoxes)]

for (x, y) in rabbitcoords:
    animals.append(Animal(x0=x, y0=y, init_energy=rabben, species=RABBIT))
for (x, y) in foxcoords:
    animals.append(Animal(x0=x, y0=y, init_energy=foxen, species=FOX))

t_rabcoordsx = []   # track the coordinates of the animals in each step in these arrays
t_rabcoordsy = []
t_foxcoordsx = []
t_foxcoordsy = []

rabbitnums, foxnums = [nrabbits], [nfoxes]  #track the number of rabbits and foxes too

animfigs = []

for i in prbar(range(steps), max_value = steps, redirect_stdout=True):          # NOTE: substitute with for i in range(steps) if progressbar2 is not installed

    # step with each animal in a random direction
    directions = np.random.randint(0, 5, size=len(animals))
    for animal, direction in zip(animals, directions):

    # generate newborn rabbits...
    rabbit_is_born_here = np.random.rand(len(coords)) <= rabbit_newborn_chance
    newrabbits = coords[rabbit_is_born_here]
    for (x, y) in newrabbits:
        animals.append(Animal(x0=x, y0=y, init_energy=rabben, species=RABBIT))

    #...  and foxes
    fox_is_born_here = np.random.rand(len(coords)) <= fox_newborn_chance
    newfoxes = coords[fox_is_born_here]
    for (x, y) in newfoxes:
        animals.append(Animal(x0=x, y0=y, init_energy=foxen, species=FOX))

    # interact if two animals are at the same coordinates
    for j, animal1 in enumerate(animals):
        for animal2 in animals[j:]:
            if (animal1.x == animal2.x and
                animal1.y == animal2.y):

    # clean up corpses
    dead_indexes = []
    for j, animal in enumerate(animals):
        if animal.isDead:
    animals = list(np.delete(animals, dead_indexes))

    # count animals and log
    foxnum, rabnum = 0,0
    for animal in animals:
        if animal.species == RABBIT:
            rabnum += 1
        elif animal.species == FOX:
            foxnum += 1
    # print(rabnum, foxnum, len(dead_indexes))

    # get and log animal coordinates
    rabcsx = []
    rabcsy = []
    foxcsx = []
    foxcsy = []
    for animal in animals:
        if animal.species == RABBIT:
            # ax.plot(, animal.y, 'bo')
        elif animal.species == FOX:
            # ax.plot(animal.x, animal.y, 'ro')


#Display the movement on an animation
fig, ax = plt.subplots()
fig.suptitle("Hunger Games")
ax.set_xlim(0, gridxsize-1)
ax.set_ylim(0, gridysize-1)

rabpc, = ax.plot(t_rabcoordsx[0], t_rabcoordsy[0], 'bo', label='rabbit')
foxpc, = ax.plot(t_foxcoordsx[0], t_foxcoordsy[0], 'ro', label='fox')

txt = ax.text(0.1, 0.1,'', ha='center', va='center', alpha=0.8,
              transform=ax.transAxes, fontdict={'color':'black', 'backgroundcolor': 'white', 'size': 10})

#initialize the animation:
def anim_init():
    rabpc.set_data(t_rabcoordsx[0], t_rabcoordsy[0])
    foxpc.set_data(t_foxcoordsx[0], t_foxcoordsy[0])
    txt.set_text('rabbits: {}\nfoxes:{}'.format(rabbitnums[0], foxnums[0]))
    return rabpc, foxpc, txt

#update the plot to the i-th frame:
def animate(i):
    rabpc.set_data(t_rabcoordsx[i], t_rabcoordsy[i])
    foxpc.set_data(t_foxcoordsx[i], t_foxcoordsy[i])
    txt.set_text('rabbits: {}\nfoxes:{}'.format(rabbitnums[i], foxnums[i]))
    return rabpc, foxpc, txt

#construct and display the animation
im_ani = animation.FuncAnimation(fig, animate, init_func=anim_init, frames=steps,
            interval=500, repeat=False, save_count=10, blit=True)

#plot population vs time
plt.plot(rabbitnums, 'b-', label="rabbits",)
plt.plot(foxnums, 'r-', label="foxes")
plt.suptitle("Population VS time")

#plot rabbuts vs foxes
plt.suptitle("Rabbit vs fox population")
plt.plot(rabbitnums, foxnums)

Sample output: grid image popt enter image description here

  • 1
    \$\begingroup\$ It would be interesting if procreation could only happen if two animals of the same species meet, costs energy, and can result in more than one offspring (with different distributions for rabbits and foxes). \$\endgroup\$ – Graipher Oct 7 '19 at 16:20
  • \$\begingroup\$ You say adding new species is fairly simple, but with the current code structure I'd have my doubts. \$\endgroup\$ – Mast Oct 8 '19 at 6:12
  • \$\begingroup\$ @Mast Right, maybe no-brainer would be a better choice of words. It indeed could be made more extendable. \$\endgroup\$ – Neinstein Oct 9 '19 at 9:49



This is suspicious. Usually you should just

#!/usr/bin/env python

and make sure that your environment uses the correct Python. Among other things, it'll make it easier for you to switch to system Python if you need it.


FOX = 1

It's good that you're assigning symbols for these values, but you should consider taking it a step further and using https://docs.python.org/3/library/enum.html .


Rather than these variables

# initial number of rabbits and foxes
nrabbits = 10
nfoxes = 5

# energy of a freshly spawned rabbit/fox
rabben = 20
foxen = 20

# chance of a new fox/rabbit being spawned at a gridpoint on a step
rabbit_newborn_chance = 0.01
fox_newborn_chance = 0.01

being tracked individually, consider making a Species class with two instances, each having the above attributes. You can then factor out a lot of common code.

Type hinting

PEP484 type hints can help, here; for instance:

def __init__(self, x0: int, y0: int, init_energy: float, species: int):

though the type for that species would change if you make it an Enum.

Variable names


should be


No-op branches

This can be deleted:

    if direction == STAY:

Global code

Your Animal class is a good start. You should move the global-scoped code starting with animals = [] into functions, as well.


# interact if two animals are at the same coordinates
for j, animal1 in enumerate(animals):
    for animal2 in animals[j:]:
        if (animal1.x == animal2.x and
            animal1.y == animal2.y):

This will get slow as you scale up to larger and larger grid sizes. Consider maintaining a nested grid list tracking where all of your animals are. List lookup is O(1) so this will be fast, and you won't have to compare all animal coordinates, which seems to be O(n^2).

Data vis

Your Rabbit vs fox population graph is cool, but it quickly becomes muddied when the system converges. This should be re-represented with another dimension, either density or time. There are many approaches to this - you could decrease the opacity of the line, or you could show a heat map, or you could make the colour of the line parametric on time and use a jet pallette so that it becomes clearer where convergence is.

Another option is to show your population vs time graph with the series smoothed and confidence intervals shown; an example is shown here - https://stackoverflow.com/questions/27164114/show-confidence-limits-and-prediction-limits-in-scatter-plot


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.