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.
imgs/parrots_alpha.pngload_image_gray()) and show it.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, :])
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.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 smallestnp.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.
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:
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_array() can take a second argument, color:
e.g. plot_array(x, color)np.tile() help you?)
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
In this diagram a,b,c,d,g,h,i are all neighbours of E
a b c
d E f
g h iThe 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
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.
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 changeblinker the 3x1 pattern above, which should switch between vertical and horizontalglider 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 simulate(cells, n):
show_image(cells)
for i in range(n):
cells = next_generation(cells)
update_image(cells)
#simulate(blinker, 10)
rake = lifeparsers.load_life("life/rake.l", padding=20)
simulate(rake, 25)
In [41]:
rake = lifeparsers.load_life("life/rake.l", padding=20)
simulate(rake, 200)
In [42]:
mystery = np.rot90(lifeparsers.load_life("life/schick.l", padding=30))
simulate(mystery, 60)
Despite the extremely simple rules, patterns in the Game of Life have an enormous diversity of behavior, with strange and intricate constructions like puffer trains, rakes, breeders, Herschel tracks and guns. It has been shown to be possible to create patterns which copy themselves (reproduce), although these are extremely complex.
It is even possible to implement a complete working computer as a pattern in the Game of Life -- a computer fabricated from cells which can compute anything a "real world" computer can (but extraordinarly slowly).
Fair warning: this is a hard problem. It is acheivable, but only attempt it if you are both enthusaistic and very comfortable with the array operations we have seen so far.
Extend the fireworks code to implement the boids animation algorithm explained at http://www.vergenet.net/~conrad/boids/pseudocode.html (a high level explanation can be read at https://cs.stanford.edu/people/eroberts/courses/soco/projects/2008-09/modeling-natural-systems/boids.html)
This models the "flocking" behaviour of birds. It is widely used in computer animation. For example, the motion of the famous wildebeest stampede scene in the Lion King was computed with a modified Boids algorithm.
Represents your boids as a single Numpy array. Use columns to represent x,y,dx,dy, and any other per-boid values you need.
Use the plot_array from the fireworks example to plot your results. I recommend calling plt.clf() before each plot_array() to clear the drawing each time (otherwise you will run into slowdowns, and you'll have a big trail behind all of your Boids).
Start with a simple 2D model of the Boids. You can implement a 3D model as well, but you will need to be able to project it onto 2D to draw it!
1D arrays can be used to represent sounds. In this task, you will implement a digital drum machine using only simple mathematical functions.
Image credit: By Eriq at Dutch Wikipedia, CC BY-SA 3.0, https://commons.wikimedia.org/w/index.php?curid=38020306
A pure tone is a sine wave. A sine wave can be computed by taking the sin() of an increasing number (e.g. generated by linspace)
In [26]:
# Warning: loud bass!
# 10000 = 10000 samples = ~1/4 second of sound (44100 samples/second)
tone = np.sin(np.linspace(0,100,10000))
play_sound(tone)
Different volume can be achieved by scaling the tone, and different pitches by scaling the value inside the np.sin()
In [27]:
# multiply *outside* the sine (try changing the 0.2)
quiet_tone = np.sin(np.linspace(0,100,10000))*0.2
play_sound(quiet_tone)
In [28]:
# multiply *inside* the sine (try changing the 4)
high_tone = np.sin(np.linspace(0,100,10000)*4)
play_sound(high_tone) # loud!
In [ ]:
# an exponentially decaying function
# this can be used to vary volume over time
a = np.linspace(0,10,200)
plt.plot(a, np.exp(-a))
You can distort sounds using np.tanh(), as shown below. Use this to improve your drum sound.
In [ ]:
tone = np.sin(np.linspace(0,100,10000))
# note we scale up the array *inside* the tanh() and then scale it down again
play_sound(np.tanh(tone*10)*0.5)
In [ ]:
# Solution goes here
In [ ]:
noise = np.random.uniform(-1,1,10000)*0.2
play_sound(noise)
Implement a hihat drum sound. Create a long and a short hihat.
In [ ]:
# Solution goes here
In [ ]:
# Solution goes here
Finally, assemble your sounds into a beat by concatenating together your drum sounds. You will need to pad your sounds with zeros such that they are all the same length if you don't want the timing to be horrible!
In [ ]:
# Solution goes here