Week 8 : Lab

Data structures: Numerical arrays

CS1P - University of Glasgow - John Williamson - 2018/2019 Python 3.x

In [1]:
from image_audio import * 
import numpy as np
%matplotlib inline
print("Import worked OK")

Import worked OK

Lab exercise

You must make a reasonble attempt at this exercise to gain a tick for this work.

Remember to save your work frequently!

Purpose of this lab

This lab will familiarise you with using numerical arrays in Python using Numpy. This includes:

  • creating arrays
  • multi-dimensional slicing
  • joining arrays
  • array arithmetic and broadcasting
  • vectorised operations
  • images as arrays

Before the lab

  • Complete at least the A exercises.

During the lab session

  • Complete the B exercise.

A: Quick problems

A.1 Box matrix

Write a function box_matrix(rows, cols) which creates a new matrix (a 2D NumPy array) with the given number of rows and columns, such that all elements are 0, except for the rows and columns at the "edges" of the matrix, which should be 1, (a "hollow" or "box" matrix").



box_matrix(2,2) =
    array([[ 1.,  1.],
           [ 1.,  1.]])

box_matrix(4,4) =
    array([[ 1.,  1.,  1.,  1.],
           [ 1.,  0.,  0.,  1.],
           [ 1.,  0.,  0.,  1.],
           [ 1.,  1.,  1.,  1.]])

box_matrix(3,6) =
    array([[ 1.,  1.,  1.,  1.,  1.,  1.],
           [ 1.,  0.,  0.,  0.,  0.,  1.],
           [ 1.,  1.,  1.,  1.,  1.,  1.]])

In [28]:
import numpy as np

def box_matrix(rows, cols):
    result = np.zeros((rows, cols))
    result[0,:] = result[0,:] + 1
    result[1:,0] = result[1:,0] + 1
    result[-1,1:] = result[-1,1:] + 1
    result[1:-1,-1] = result[1:-1,-1] + 1
    return result
print(box_matrix(4, 4))

[[1. 1. 1. 1.]
 [1. 0. 0. 1.]
 [1. 0. 0. 1.]
 [1. 1. 1. 1.]]

A. 2 Tiling

np.tile(a, shape) will copy a multiple times in a grid shape. For example, np.tile(a, (2,2)) will copy an array twice in columns and twice in rows, joining the result together:

In [29]:
a = np.array([[1,0], [0,4]])
np.tile(a, (2,2))

array([[1, 0, 1, 0],
       [0, 4, 0, 4],
       [1, 0, 1, 0],
       [0, 4, 0, 4]])

Using np.tile create an array of the numbers [1,2,3,4] tiled into 4 columns and 4 rows, so that the resulting output is a 4x16 array of values, like this:


Show the result as an image using show_image(x), but divide every element of x by 4 before showing it.


In [38]:
a = np.array([1, 2, 3, 4])
show_image(np.tile(a, (4,4)) / 4.0)

A.3 Alpha mask

An grayscale image consists of a 2D array (matrix) of brightness values. Combining two images together can be done by adding them together and dividing by 2.0, as below.

In [39]:
parrots = load_image_gray("imgs/parrots.png")
plane = load_image_gray("imgs/plane.png")

These are just plain old arrays, and we can do anything we want to the elements

In [40]:
show_image(plane*2) # increase contrast
show_image(np.abs(plane - parrots)) # show absolute difference of images

# the equivalent of Photoshop's "Darken"  mode
show_image(np.where(plane<parrots, plane, parrots))

But to composite on image on top of another (e.g. to place an actor filmed in front of a green screen into a CGI background), we need to selectively mix to gether images. This is called alpha masking (this why the transparency option in some programs is called alpha).

There is a simple process: you take two images $X$ and $Y$, and a special third image $A$ called the alpha mask or sometimes the matte. $A$ is a grayscale image; where $A$ is white (elements = 1.0), $X$ shows through. Where $A$ is black (elements = 0.0), $Y$ shows through. Where $A=0.5$, $0.5X+0.5Y$ is mixed together.

The general formula for this masking is very simple $$Z = AX + (1-A)Y,$$ where $Z$ is the output image.

  • An alpha mask is contained in the file imgs/parrots_alpha.png
  • Load it (use load_image_gray()) and show it.
  • Compute the alpha composited image (parrots on background) by implementing the equation above using ndarrays. This is very simple and should not involve a loop!
  • Show the composited image.

Optional extension

Make this work in colour (use load_image_colour() to load a colour image; show_image() will still work to show it)

In [55]:
# Z = AX + (1-A)Y
alpha = load_image_gray("imgs/parrots_alpha.png")
parrots = load_image_gray("imgs/parrots.png")
beach = load_image_gray("imgs/beach.png")
result = alpha * parrots + (1 - alpha) * beach

A.4 Sunspots

Astronomers have counted the number of visible sunpots on the sun every day since 1749. The code below loads this data (average daily sunspot count over one month intervals) as a Numpy array. It has three columns:

index, year, mean sunspots

year is in the date in fractional years (e.g. 1800.5 means half way through 1800).

You can think of this like an Excel spreadsheet with data in it.

The index just increments from 1 upwards, and isn't interesting to us.

In [56]:
# load the data
sunspots = np.loadtxt("sunspots.csv", delimiter=",")    
print((sunspots.shape, sunspots.dtype))

((2820, 3), dtype('float64'))

In [57]:
# print the first five rows of this data

print(sunspots[:5, :])

[[   1.         1749.           58.        ]
 [   2.         1749.08333333   62.6       ]
 [   3.         1749.16666667   70.        ]
 [   4.         1749.25         55.7       ]
 [   5.         1749.33333333   85.        ]]

We can plot the sunspot data directly with plot. We will see the characteristic 11 year solar activity cycle.

In [58]:
# plot second column (time in years) against third (number of sunspots)
plt.plot(sunspots[:,1], sunspots[:,2])
plt.ylabel("Avg. daily sunpots")

Text(0, 0.5, 'Avg. daily sunpots')

There are two useful operations in Numpy:

  • np.argmin(x) returns the index at which x is smallest
  • np.argmax(x) returns the index at which x is largest

In [59]:
x = [0.5, -0.2, 0.9, 1.8]


In [60]:
# we can use the numbers returned as indices


Without using any loops, answer the following questions using array operations. You should be able to answer these in a single line of code.

  1. What is the most number of sunspots seen?
  2. What is the least number of sunspots seen?
  3. What is the average (mean) number of sunspots seen each month?
  4. What is the average number of sunspots for March? Hint: there are 12 months in the year!
  5. In which year were most sunspots seen?

In [99]:
import math

print("1. Most sunspots: ", sunspots[np.argmax(sunspots[:,-1]), -1])
print("2. Least sunspots: ", sunspots[np.argmin(sunspots[:,-1]), -1])
print("3. Average sunspots: ", round(np.mean(sunspots[:,-1]), 2))
print("4. Average in March: ", round(np.mean(sunspots[2::12, -1]), 2))
print("5. Most sunspots seen in: ", math.floor(sunspots[np.argmax(sunspots[:,-1]), 1]))

1. Most sunspots:  253.8
2. Least sunspots:  0.0
3. Average sunspots:  51.27
4. Average in March:  50.04
5. Most sunspots seen in:  1957

B: Simulations with arrays

You have been hired to demonstrate the artistic and expressive aspects of computer science. You have to build an installation to be displayed in a museum to show off what can be done.

You have two animated exhibits to produce:

  • Firework particle explosion
  • The Game of Life

You will implement this using numerical arrays to perform computations quickly and elegantly.

B.1 Firework displays

We missed Guy Fawkes day, but we can simulate fireworks using a simple particle system. We will use very (very) basic physics to simulate how particles (in this case, little bits of fireworks) move around. We will use numerical arrays (NumPy arrays) to compute the movement of a "whole" firework at once.

The output from your code for this part should look like this:

We'll simulate this using a very set of simple equations that apply to a whole set of particles at the same time (i.e. vectorized operations).

Particle system

Assume we can model a single firework explosion at a given time as a set of N points in 2D, i.e. an array of N rows of 2 columns (x and y).

[[x, y]
 [x, y]
 [x, y]

To simulate a firework exploding, we need to simulate the motion of these points "bursting out" from a central position and falling down.

Assume we also have some information about the initial velocity of each particle (i.e. which way it went at the moment of detonation). This is a vector, and has two components, "dx" and "dy", the rate of movement in each of those two directions.

We can write a 2D vector in space at time $t$ as $\textbf{x_t}=[x,y]$ and a 2D velocity vector as $\textbf{dx_t} = [dx, dy]$.

We can predict where the particle is at time $t$ from where it was in the last time step $t-1$ using the equations: $$\textbf{x}_{t} = \textbf{x}_{t_1} + \textbf{dx}_{t_1}$$ $$\textbf{dx}_{t} = d \textbf{dx}_{t-1} - \textbf{g}$$

$d$ is a damping factor (from air resistance). $g$ is the gravitational vector, which accelerates towards the ground.

Note that x, g and dx are vectors (they have $x$ and $y$ components), not just plain numbers!

Random numbers

We can use random numbers to give initial values for the velocities (after all, the particles all spread out in random ways when a firework goes off). The function np.random.normal(mean, std, shape) generates normally distributed values, centered on mean with a "spread" given by std, in an array of the given shape.

In [100]:
# 4 element vector, centered on 0, spread of 0.1
print((np.random.normal(loc=0, scale=0.1, size=(4,))))

[ 0.02838059 -0.2023798  -0.14933385  0.14830499]

In [101]:
# 8 element vector, centered on 0, spread of 10
print((np.random.normal(loc=0, scale=10, size=(8,))))

[ -3.38668139   6.03132396 -10.45926129 -13.24456308  -6.43303956
 -12.26632799  -5.64600167  13.98906436]

In [102]:
# 6x6 matrix, centered on 0, spread of 2
print((np.random.normal(loc=0, scale=2, size=(6,6))))

[[-1.46296018 -0.11015676  3.35095462  0.57312552  2.36993829 -0.08212527]
 [-2.32524807 -1.64108725  3.77894467  0.53982246 -0.06367447  2.16395739]
 [-1.20792648  3.05398352 -0.70815021  0.54563987  1.84075902  0.33854205]
 [ 0.25611395  1.57044427  0.00949805 -0.82272334 -1.53736211 -0.24145984]
 [ 4.37385311  0.58552502  1.62515387 -1.88049632  1.19627449  0.07664221]
 [-1.06445536 -2.85427237  0.51187776 -0.58972294 -1.70010902 -1.95125376]]

Task 1

Implement the equations of motion for a single firework.

  • You should start with all points at the same position $\textbf{x}$ (e.g. all at (0,0)), but with random velocities, using np.random.normal(), as above.

  • Start with 50 particles in the firework. Remember, the positions of every particle should be in one single array!

  • To implement the gravity, remember you can modify a slice of an array (e.g. a column slice).

  • Run the equations for 50 steps, updating the position and velocity. Each iteration should compute a new array.

  • Use plot_array(pts) defined below to plot the positions on each timestep. Use the s optional parameter to set the size of the points to 50-i, where i is the iteration number.

  • To start, use gravity g=[0, -0.01], damping d=0.9. Adjust these as you feel best.

  • At the end, you should get a single, monochrome firework like this:

In [126]:
# 𝐱𝑡 = 𝐱𝑡1 + 𝐝𝐱𝑡1
# 𝐝𝐱𝑡 = 𝑑𝐝𝐱𝑡-1 − 𝐠

def plot_array(pts, c=None, s=8):
    """Plot pts (must be a Nx2 array) as points. 
    If c is specified, it must be an Nx4 array of RGBA values, giving
    the colour for each point."""
    if c is None:
        plt.scatter(pts[:,0], pts[:,1], c='r', edgecolor="none", s=s)
        plt.scatter(pts[:,0], pts[:,1], c=np.clip(c,0,1), edgecolor="none", s=s)

In [128]:
g = [0, -0.01]
d = 0.9
velocities = np.random.normal(0,1,(50,2))
points = np.zeros((50,2))

for i in range(50):
    plot_array(points, s=50-i)
    points += (velocities*d)
    velocities = (velocities*d) + g

KeyboardInterrupt                         Traceback (most recent call last)
<ipython-input-128-83f7d5798635> in <module>
      6 print(points)
      7 for i in range(50):
----> 8     plot_array(points, s=50-i)
      9     points += (velocities*d)
     10     velocities = (velocities*d) + g

<ipython-input-126-4382a5ff709c> in plot_array(pts, c, s)
      8     plt.gca().set_facecolor('black')
      9     if c is None:
---> 10         plt.scatter(pts[:,0], pts[:,1], c='r', edgecolor="none", s=s)
     11     else:
     12         plt.scatter(pts[:,0], pts[:,1], c=np.clip(c,0,1), edgecolor="none", s=s)

/usr/local/lib/python3.7/site-packages/matplotlib/pyplot.py in scatter(x, y, s, c, marker, cmap, norm, vmin, vmax, alpha, linewidths, verts, edgecolors, plotnonfinite, data, **kwargs)
   2845         verts=verts, edgecolors=edgecolors,
   2846         plotnonfinite=plotnonfinite, **({"data": data} if data is not
-> 2847         None else {}), **kwargs)
   2848     sci(__ret)
   2849     return __ret

/usr/local/lib/python3.7/site-packages/matplotlib/__init__.py in inner(ax, data, *args, **kwargs)
   1599     def inner(ax, *args, data=None, **kwargs):
   1600         if data is None:
-> 1601             return func(ax, *map(sanitize_sequence, args), **kwargs)
   1603         bound = new_sig.bind(ax, *args, **kwargs)

/usr/local/lib/python3.7/site-packages/matplotlib/axes/_axes.py in scatter(self, x, y, s, c, marker, cmap, norm, vmin, vmax, alpha, linewidths, verts, edgecolors, plotnonfinite, **kwargs)
   4525                 self.set_ymargin(0.05)
-> 4527         self.add_collection(collection)
   4528         self.autoscale_view()

/usr/local/lib/python3.7/site-packages/matplotlib/axes/_base.py in add_collection(self, collection, autolim)
   1872         if autolim:
-> 1873             self.update_datalim(collection.get_datalim(self.transData))
   1875         self.stale = True

/usr/local/lib/python3.7/site-packages/matplotlib/collections.py in get_datalim(self, transData)
    200                 transform.frozen(), paths, self.get_transforms(),
    201                 offsets, transOffset.frozen())
--> 202             result = result.inverse_transformed(transData)
    203         else:
    204             result = transforms.Bbox.null()

/usr/local/lib/python3.7/site-packages/matplotlib/transforms.py in inverse_transformed(self, transform)
    520         of *transform*.
    521         """
--> 522         return self.transformed(transform.inverted())
    524     coefs = {'C':  (0.5, 0.5),

/usr/local/lib/python3.7/site-packages/matplotlib/transforms.py in inverted(self)
   2403     def inverted(self):
   2404         # docstring inherited
-> 2405         return CompositeGenericTransform(self._b.inverted(), self._a.inverted())

/usr/local/lib/python3.7/site-packages/matplotlib/transforms.py in inverted(self)
   2403     def inverted(self):
   2404         # docstring inherited
-> 2405         return CompositeGenericTransform(self._b.inverted(), self._a.inverted())

/usr/local/lib/python3.7/site-packages/matplotlib/transforms.py in inverted(self)
   1807             if self._shorthand_name:
   1808                 shorthand_name = '(%s)-1' % self._shorthand_name
-> 1809             self._inverted = Affine2D(inv(mtx), shorthand_name=shorthand_name)
   1810             self._invalid = 0
   1811         return self._inverted

<__array_function__ internals> in inv(*args, **kwargs)

/usr/local/lib/python3.7/site-packages/numpy/linalg/linalg.py in inv(a)
    549     signature = 'D->D' if isComplexType(t) else 'd->d'
    550     extobj = get_linalg_error_extobj(_raise_linalgerror_singular)
--> 551     ainv = _umath_linalg.inv(a, signature=signature, extobj=extobj)
    552     return wrap(ainv.astype(result_t, copy=False))


Task 2

  • Plot several fireworks, each with a different random starting position.
  • Start with 5 fireworks.
  • plot_array() can take a second argument, color: e.g. plot_array(x, color)
  • The colours given to plot_array must be an Nx4 array (red, green, blue, alpha). The alpha column should initially be set to 1, so that the fireworks don't start out transparent. RGB should be values between 0 and 1.
  • Choose a random colour for each firework, but with the alpha component set to 1.0. Note that you need to give a colour for each particle in the firework. (hint: how can np.tile() help you?)
  • Make the colours fade as the simulation runs (i.e. they should get darker or more transparent in future time steps). Think about how the damping changed the velocity.

In [229]:
def plot_firework(particles, start_x=0, start_y=0):
    g = [0, -0.01]
    d = 0.9
    # Generate firework colour
    colour = np.tile(np.random.randint(low=0, high=100, size=4) / 100.0, (50, 1))
    a = 1.0
    for i in range(particles):
        colour[i,-1] = a
        a *= d
    velocities = np.random.normal(0,1,(50,2))
    points = np.zeros((50,2))
    points[:,0] = start_x
    points[:,1] = start_y
    for i in range(50):
        plot_array(points, colour, s=50-i)
        points += (velocities*d)
        velocities = (velocities*d) + g
plot_firework(50, 2, 2)

B.2 Game of Life

The Game of Life is a very famous simulation; a kind of cellular automaton. It is not a game in the traditional sense of the word, but a kind of simulation which has many interesting and surprising properties. Given a particular starting patterns, unexpectedly complicated behaviours evolve, including moving objects ("spaceships"), repeating patterns ("oscillators"), patterns that shoot out spaceships ("guns") and many other structures.

Oscillating patterns

A large, slow spaceship

A "glider gun", shooting out little glider spaceships

A "fuse reaction"

The Game of Life works by simulating the birth and death of a very simplified model of a 2D grid of "cells" according to really simple rules. All of the behaviours above arise from these rules. It was invented by the mathematician John H. Conway in 1969

John H. Conway [Image Credit "Thane Plambeck" - "http://www.flickr.com/photos/thane/20366806/", CC BY 2.0, https://commons.wikimedia.org/w/index.php?curid=13076802]

The Game of Life assumes that these "cells" are represented as a big 2D array (like a giant chessboard), where each element of the array can be 1 (on/alive) or 0 (off/dead). Rules are defined which, given a array of cells in one instant, determine how the cells should evolve in the next instant.

There's more background in this blog post


  • We imagine that we loop over every entry in an array, and at that position look at the 3x3 neighbourhood of cells that surrounds that entry. That is every array entry that is one step horizontally, vertically or diagonally away from the point we are considering.

In this diagram a,b,c,d,g,h,i are all neighbours of E

a b c
d E f
g h i
  • For each cell in the array:
    • Check all eight neighbours in a 3x3 grid around the cell
    • If cell was already on, and there are exactly 2 neighbours, the cell stays on in the next timestep.
    • If there are exactly three neighbours, the cell becomes on in the next timestep, regardless of whether it was on or off before.
    • Otherwise, the cell becomes off in the next timestep.

The eight neighbours of a cell look like this:

The rule

Still Lifes

The Game of Life rule has very interesting properties. Some arrangements of cells ("patterns") are stable and never change, like 4 on cells in a 2x2 grid, surrounded by all off cells.

0 0 0 0
0 1 1 0
0 1 1 0
0 0 0 0

These are called still lifes.


Other patterns continuously oscillate between two states, like 3 cells in a line:

0 0 0 0 0
0 1 1 1 0
0 0 0 0 0

Patterns like these are called oscillators


Remarkably, there are discrete equivalents of waves in this simulation; patterns which appear to move across empty space.

   0 0 0 0 0
   0 1 1 1 0
   0 1 0 0 0
   0 0 1 0 0
   0 0 0 0 0

These are called spaceships and are much, much rarer than any of the other patterns we have seen so far.


The Game of Life is very easy to represent as a big numerical array, each element of the array corresponding to a cell. Your task is to implement the Game of Life using arrays.


You need to write a function next_generation(cells), which given a numerical array cells consisting of 0s and 1s, returns a new array of 0s and 1s obeying the Game of Life rules.

You can use show_image(cells) to show the output array. Test your function on the block and blinker examples before going further.

Then, write a function simulate(cells, n) which will apply next_generation n times to the cells, and show an animation of the process by redrawing the image on each generation.

You can use update_image(cells) to update a drawing without creating a completely new image for each update; this lets you show an animation as simulate() runs. Note that you must call show_image() the first time you draw the image (to create a new blank figure canvas), and call update_image() to update the figure.

Correct simulation

Note: to get the simulation correct, you need create a new array and copy the values in. If you try and modify the cells in place, you will get completely incorrect behaviour (unless you are very clever with the array ops!).


Keep the shape of cells the same on each application of next_generation(). Note that you will have to do something at the edges of the array, where there aren't a full set of neighbour cells (e.g. what is top left of the the top left cell?).

Wrapping at the edges

A common solution is just to assume that the array "wraps around" at the edges (e.g. that the leftmost element is next to the rightmost element). This is very easy to implement using the % operator, which allows you to force a value to lie within a certain range; see the example below to see why:

In [230]:
# forwards:
for i in range(12):
    print(i, i % 5, " ", end=' ')

# backwards
for i in range(0,-12,-1):
    print(i, i % 5, " ", end=' ')

0 0   1 1   2 2   3 3   4 4   5 0   6 1   7 2   8 3   9 4   10 0   11 1   
0 0   -1 4   -2 3   -3 2   -4 1   -5 0   -6 4   -7 3   -8 2   -9 1   -10 0   -11 4   

Test patterns

There are some test patterns provided for you. load_pattern(name, padding) returns a 2D array with the given pattern, with padding 0 (off) cells around it. You can check that this is just a standard Numpy array. The padding argument allows you to get an array with enough room for a pattern to "grow".

These patterns include:

  • block the 2x2 pattern above, which should not change
  • blinker the 3x1 pattern above, which should switch between vertical and horizontal
  • glider the spaceship above, which should move up and to the right.

There are many more.

In [239]:
import lifeparsers

# load the block
block = lifeparsers.load_life("life/block.l")

# load the blinker
blinker = lifeparsers.load_life("life/blinker.l")

# load the glider (with some extra padding)
glider = lifeparsers.load_life("life/glider.l", padding=4)

In [40]:
import lifeparsers

# load the blinker
blinker = lifeparsers.load_life("life/blinker.l")

#Check all eight neighbours in a 3x3 grid around the cell

#If cell was already on, and there are exactly 2 neighbours, the cell stays on in the next timestep.

#If there are exactly three neighbours, the cell becomes on in the next timestep,
#regardless of whether it was on or off before.

#Otherwise, the cell becomes off in the next timestep.

def next_generation(cells):
    rows = len(cells[:, 0])
    cols = len(cells[0,:])

    new_cells = np.zeros((rows, cols))
    for y in range(rows):
        for x in range(cols):
            neighbours = [[(x-1)%cols, y], [(x+1)%cols, y], [x, (y+1)%rows], [x, (y-1)%rows],
                          [(x+1)%cols, (y+1)%rows], [(x-1)%cols, (y-1)%rows],
                          [(x+1)%cols, (y-1)%rows], [(x-1)%cols, (y+1)%rows]]
            if sum(cells[y, x] for x, y in neighbours) == 2 and cells[y, x] == 1:
                new_cells[y, x] = 1
            elif sum(cells[y, x] for x, y in neighbours) == 3:
                new_cells[y, x] = 1
                new_cells[y, x] = 0
    return new_cells

def simulat