```
In [1]:
```## RUN ME FIRST!
from image_audio import *
import numpy as np
np.set_printoptions(suppress=True)
%matplotlib inline
print("Import worked OK")

```
```

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").

**DO NOT USE ANY FORM OF LOOP**

So

```
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))

```
```

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

```
Out[29]:
```

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:

```
[[1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4],
[1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4],
[1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4],
[1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4]]
```

Show the result as an image using `show_image(x)`

, but divide every element of `x`

by 4 before showing it.

DO NOT USE ANY FORM OF LOOP.

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

```
```

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

```
```

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.

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
show_image(result)

```
```

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))

```
```

```
In [57]:
```# print the first five rows of this data
print(sunspots[:5, :])

```
```

`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.xlabel("Year")
plt.ylabel("Avg. daily sunpots")

```
Out[58]:
```

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]
print(np.argmin(x))
print(np.argmax(x))

```
```

```
In [60]:
```# we can use the numbers returned as indices
print(x[np.argmin(x)])
print(x[np.argmax(x)])

```
```

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

- What is the most number of sunspots seen?
- What is the least number of sunspots seen?
- What is the average (mean) number of sunspots seen each month?
- What is the average number of sunspots for March?
*Hint: there are 12 months in the year!* - 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]))

```
```

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.

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).

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!

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,))))

```
```

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

```
```

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

```
```

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."""
plt.gca().set_facecolor('black')
if c is None:
plt.scatter(pts[:,0], pts[:,1], c='r', edgecolor="none", s=s)
else:
plt.scatter(pts[:,0], pts[:,1], c=np.clip(c,0,1), edgecolor="none", s=s)
plt.xlim(-10,10)
plt.ylim(-10,10)
display.display(plt.gcf())
display.clear_output(wait=True)
time.sleep(0.05)

```
```

```
In [128]:
```g = [0, -0.01]
d = 0.9
velocities = np.random.normal(0,1,(50,2))
points = np.zeros((50,2))
print(points)
for i in range(50):
plot_array(points, s=50-i)
points += (velocities*d)
velocities = (velocities*d) + g

```
```

- 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)
plot_firework(50, 2, 2)

```
```

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 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
```

**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.

`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.

**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?).

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=' ')
print()
# backwards
for i in range(0,-12,-1):
print(i, i % 5, " ", end=' ')

```
```

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")
show_image(block)
# load the blinker
blinker = lifeparsers.load_life("life/blinker.l")
show_image(blinker)
# load the glider (with some extra padding)
glider = lifeparsers.load_life("life/glider.l", padding=4)
show_image(glider)

```
```

```
In [40]:
```import lifeparsers
# load the blinker
blinker = lifeparsers.load_life("life/blinker.l")
#show_image(blinker)
#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
else:
new_cells[y, x] = 0
return new_cells
def simulat