PX1224: Computational Skills for Problem Solving

Two Dimensional Arrays

Introduction

This week is concerned with 2 dimensional (2-d) arrays. 2-d arrays are useful as they allow more complex data to be represented in the computer, such as matrices and images. They also allow related data (eg extension and force in a Hooke's law experiment) to be kept aligned. In this week's work, you will create some 2-d array and matrices and do some calculations, reinforcing material that you have already met in the maths classes.

Syntax Summary

Defining 2-d Arrays

Function Syntax
Define a 2-d array a = array([[1,2],[3,4]])
Define a 2-d array of zeros z = zeros( (10, 10) )
Define a 2-d array of ones o = ones( (10, 10) )
Define a 2-d grid of numbers y, x = mgrid[0:5,0:5]
Stack 1-d arrays to make a 2-d array a = column_stack([x,y])

Defining 2-d Arrays

Function Syntax
Return an entry of a 2-d array a[2,2]
Return part of a 2-d array a[0:2, 0:2]
return entries with a step greater than one a[0:6:2, 0:4:2]
1d slice of a 2-d array a[:,0]
a[0,:]
return all entries as a 1-d array a.flatten()

Matrices

Function Syntax
generate a matrix a = matrix([[1,2],[2,3]])
invert a matrix a.I
tranpose a matrix a.T
eigenvalues and eigenvectors of a matrix values, vectors = eig(a)

Plotting 2-d arrays

Function Syntax
make a contour plot of 2-d data contourf(x,y,z,100)
display a 2-d array imshow(a)
Show the colour bar colorbar()

Worksheet

Two dimensional arrays

A one-dimensional array is an indexed "list" of numbers, so each "block" in the array holds a number. A two-dimensional array extends this concept so that each "block" in a one-dimensional array holds another one-dimensional array.

A 1-d array:

Index 0 1 2 3 4 5 6 7 8
Value 12.3 2.3 5.2 77.3 11.2 55.2 74. 22. 99.

To input this into Python we would type

y = array(  [12.3, 2.3, 5.2, 77.3, 11.2, 55.2, 74., 22., 99]  )

Note:Above, the input to the array() function is a list of numbers (denoted by square brackets)

As a 2-d array is "a list of arrays", where each row is an array and it is well represented on a 2-dimensional table

0 1 2
0 12.3 2.3 5.2
1 22.3 23.2 86
2 -2 .3 33

Generating 2-d arrays

As with 1-d arrays, there are several ways to generate a 2-d array:

  1. Explicitly fill the elements you desire. For example, to generate the 2-d array given above in python, you would type
    d2 = array([
                  [12.3, 2.3, 5.2],
                  [22.3, 23.2, 86.],
                  [-2, .3, 33]
                ])
  2. We can also provide a 2D array of zeros or ones for any given number of rows and columns (eg, 32 rows and 40 columns) using NumPy's zeoros() and oness() functions:
    z = zeros( (32, 40) )
    o = ones( (32,40) )
  3. The function mgrid() can be used to make 2-d arrays in a similar way that arange() can be used for 1-d arrays. mgrid() will return two 2-d arrays, one that contains the x-coordinates and the other the y-coordinates. For example:
    y,x = mgrid[0:5, 0:5]
    print(x)
    print(y)
    Note: we have deliberately asked for y,x = mgrid[0:5, 0:5] and not x,y = ... Be careful with this.
    This will give the x array will have 5 rows each containing [0 1 2 3 4], while the y array will have a row of [0 0 0 0 0] followed by [1 1 1 1 1], etc. This is like a two-dimensional version of arange where arange(0, 5) would give [0 1 2 3 4].
    As with arange, you can vary the step size (the default is 1), so
    y,x = mgrid[0:5:0.1, 0:5:0.1]
    print(x)
    print(y)
    will return two arrays that have 50x50 entries. x will have values, equally spaced between 0 and 5, in steps of 0.1, in the direction of rows. Similarly, y will have values between 0 and 5 but in the direction of columns. In the examples here, we have chosen to have a square array (same number of points in x and y), but that is not necessary.
  4. mgrid() can also give arrays with a specified number of points, similar to linspace(). For example, to generate an array with 10 points in the x direction and 6 in the y direction:
    y,x = mgrid[0:5:10j, 0:5:6j]
    print(x)
    print(y)

Accessing entries in 2-d arrays

Similarly to the 1D arrays, 2D array elements can be addressed individually, or subsections can be addressed, both to set values or return values. Here, we introduce the ways of interrogating a 2-d array, many of which should be familiar from the methods we learned for 1-d arrays previously.

  1. You can print the contents of a 2d array using print
    print(d)
    
  2. There are several different ways to understand the size of the array. The size command will give you the total number of entries in the array:
    size(d)
    
    The shape command tells you the length of each of the dimensions in the array, for example
    shape(d)
    
    The len command gives the length of the first dimension in the array -- it gives no information about any additional dimensions
    len(d)
    
  3. As with 1d arrays, you can access any entry you like, The order should be familiar from matrices where M21 identifies the entry in the second row, first column. The difference with python is that the numbering starts at zero. For example
    print d2[0,0] # Value in first row, first column
    print(d2[2,0]) # Value in third row, first column
    print(d2[0,2]) # Value in first row, third column
    
  4. You can also access entries starting from the end, e.g.
    print(d2[-1,-1]) # Value in last row, last column
    print(d2[-1,0]) # Value in last row, first column
    
  5. You can access a range of entries by specifying the range that you are interested in, for example
    b2 = d2[0:2,0:2] #Values from the first 2 rows and first 2 columns (2x2)
    print(b2)
    
  6. More generally, you can also specify every other entry, every 5th entry, etc. In the 1D case, you used the syntax [start:stop:step]. In the 2D case, the syntax is [start_row:stop_row:step_row, start_column:stop_column:step_column]. As an example,
    c = d2[ : :2 , : :2]
    print(c)
    Note:In the above example only the step size of 2 is specified - the start and the end rows/columns are not specified, in which case they are assumed to be the first and last column of the 2d array.

Slicing Arrays

A 2D array can be "sliced" to return a 1D array containing one row or column of the original. So, for example, the zeroth (Python starts array with index 0) row or column is given by

r0 = d2[0,:] #Extract the first Row and call it r0
c0 = d2[:,0] #Extract the first Column and call it c0
c3 = d2[:,3] #Extract the 4th Column and call it c4

It is often useful to extract 1-dimensional slices of 2-d arrays for further investigation.

Calculations with 2d arrays

2-d Arrays can be multiplied just as 1-d arrays can be; all mathematical operations that work with single variables (i.e. single numbers) (addition, subtraction, division, raising to a power, sin etc.) work with 2D arrays element-by-element.

It is possible to find the maximum entry in a 2-d array, however, if you try using python's max() function, you will get an error

max(d2) #Will raise Error

ValueError: The truth value of an array with more than one element is ambiguous. Use a.any() or a.all(). This is because Python's max() function cannot work with NumPy 2d array objects.

To calculate the max value in a 2d NumPy array, you have to instead use the array's max() method:

d2.max() # This will work

You can sum the entries of a 2-d array using sum, and also sum up just the rows or columns by specifying the axis you want to sum along

sum(d2) # add all entries of d2
sum(d2, axis = 0) # sum entries in each column separately
sum(d2, axis = 1) # sum entries in each row separately

Truth testing also works with 2-d arrays, for example

print(d2 > 5)
print(d2[d2 > 5])

Note that the first line returns a 2-d array of True/False while the second line returns a 1-d array of values. You can also cut one array based on the values contained in a second array, for example

y,x = mgrid[0:5, 0:5]
x[y==0] #Extract those x values, when the corresponding y values are equal to 0

The above code snippet is further described by the diagram below

Reading and Writing 2-d arrays

Reading a 2-d array from a file

In many instances, data in a file or spreadsheet will have several columns, and this is more likely to be the kind of data that you generate in a lab. Last week, we learned how to read and write 1-d data, here we look at 2-d data.

As an example, we have made a file with the data that was used previously for the Hooke's law graph. It's in the file "Hooke's law data.xls", and again saved as both tab (.txt) and comma separated (.csv) files. Opening up any of these in Excel, they all look the same. But now, if you look in Notepad, you'll see some differences. In the .txt file, the columns are separated by tabs while the .csv file has comma separations.

Either tab or comma separated files can be read using loadtxt(), you just have specify the delimiter (comma, tab, space, other) so that the computer knows how to break up the lines.

spring = loadtxt( "Hookes Law Data.txt", skiprows=1)
print(spring)
spring2 = loadtxt( "Hookes Law Data.csv", skiprows=1, delimiter=',' )
print(spring2)

Note: In the first case, the delimiter isn't needed as a tab (or any whitespace) is the expected format.

Looking at spring in the variable explorer, you'll see that it's a 2-d array with 7 rows and 2 columns.

Alternatively, we can use one more trick to make loadtxt() return a set of 1-d arrays. To do this, we use the "unpack=True" parameter as

mass, length = loadtxt( "Hookes Law Data.txt", skiprows=1, unpack=True )
print(mass)
print(length)

This has unpacked each column into a separate 1-d array.

The loadtxt() function has other useful options, such as setting the data type or only using certain columns. We will introduce these in later weeks.

Writing several columns of data

You can use savetxt() to write a 2-d array to file. To write the spring data to a new file, simply do the following

savetxt("new_hookes.txt", spring, "%.2f")

If you have "unpacked" the data to 1-d arrays when reading it, you need to "re-pack" (i.e stack the 2 1D arrays in to a 2D matrix) it when writing. This is done using the column_stack() function

column_stack([mass,length])

So, to put everything together

savetxt("new_hookes.txt", column_stack([mass,length]), fmt="%.2f\t%.2f" )

The fmt option specifies the format of each line to be a float printed to two decimal places followed by a tab and a second float printed to two decimal places.

Plotting 2-d data

A 2-d array can represent a "bitmap" image. Each element in the 2d array can represent a brightness, it is easy to imagine creating a grayscale image in this way. Full colour images actually store 3 values, "redness", "blueness" and "greeness" for each value, so can be thought of as three 2D arrays. Alternatively, a 3D array, eg 1200 elements by 800 elements (ie a screen resolution of 1200x800) by 3 elements can represent a full colour image on the screen.

To generate an image like this, we first want to create x and y arrays using mgrid(), then calculate a "height" value, z, at each point in x and y. Finally, we plot a contour plot using the values of z to give the colours. The code to do this is given below

Python can show such images. The example code below produces a large 2D array, and then produces a pattern in it.

y,x = mgrid[-1: 1: .01  , -1 : 1: .01]
z = 3*y*(3*x**2-y**2)/4 + .5*cos(6*pi * sqrt(x**2 +y**2) + arctan2(x,y))
contourf(x,y,z,100)
colorbar()

Figure 1: Example of a filled contour plot, generated using the code above.

Note: The value 100 in the contourf() function determines the number of different coloured contours to plot. If you reduce this to, for example, 10 then the plot becomes less smooth.

Note: As an alternative to contourf(), you can use imshow().

imshow(z, extent=[-1,1,-1,1], origin='lower', interpolation='nearest' )

Here, you need to specify the ranges of x and y using extent, and also ask the origin to be in the lower left (default is upper left). We have also turned off interpolation of data: this plots the data as it is in the array, rather than smoothing it first.

Creating surface plots

Another way to display 2-d arrays is as 3-dimensional surface plots. Matplotlib, the plotting library we use in python, is capable of creating a variety of 3d plot types, examples of which can be seen here. One such example is shown in the following code.

from mpl_toolkits.mplot3d import Axes3D
#Create some data to plot. 
y,x = mgrid[-2.:2.:.1, -2.:2.:.1]
z = exp(-x**2 - y**2)
 
# plot data (x,y,z)
fig = plt.figure()
ax = fig.gca(projection='3d')
surf = ax.plot_surface(x, y, z, rstride=2, cstride=2, shade=False, cmap="jet")
draw()
surf.set_edgecolors(surf.to_rgba(surf._A))
surf.set_facecolors("white")
ax.set_xlabel('X Label')
ax.set_ylabel('Y Label')
show()

Figure 2: Example of a coloured wireframe surface plot, generated using the code above.

Note: use parameters rstride=n, cstride=n to plot every nth element of the 2d-array. This is handy when you try to make surface plots of very large arrays.

Note: The 3D plots are interactive. You can rotate the view by click and dragging on the plot.

Matrices

NumPy knows about a special type of 2D array, the matrix. A matrix can be set up in the same way as a 2D array, just with the function matrix() rather than array(). An alternative method also exists, which requires a little less typing: So, a 2D matrix can be generated with

#Using lists
a = matrix([  [1,2]  ,   [2,3]   ])

#Or using a string of space separated values, semicolon separated rows
b = matrix('1   -2  ; 2   3 ')

Once a matrix has been entered, NumPy allows you to multiply matrices, multiply vectors (ie a 1D array!) by a matrix, determine the transpose, determinants and inverse. We will introduce all of these in turn below.

As with all arrays in Python, mathematical operations can be performed on 2D arrays and matrices, but only if they are of correct shapes. (Shape is the python name for the dimensions of the array.) Attempting to combine (addition, subtraction etc) two arrays that are not compatable results in the error message ValueError: shape mismatch: objects cannot be broadcast to a single shape.

Matrix Addition

We can add two matrices as

a + b

Matrix Multiplication

We can multiply a matrix by a constant using

c = 5 * a
print(c)

and multiply matrices as

a * b

Try this with both arrays and matrices, and make sure you understand the difference. To do this with arrays, you could either re-enter the matrices above as arrays, or simply do

array(a) * array(b)

Inverse

We can calculate the inverse of the matrix using

print(a**-1)

and check that a multiplied by its inverse gives the identity matrix

print(a * a**-1)

Determinant

We can compute the determinant using

det(a)

Transpose

Determine the transpose

a.T

Note: The transpose is a method of the class of matrices. Recall that you met methods last week. There are a number of useful methods in the matrix class

Method Result
a.T Transpose
a.I Inverse
a.H (complex) conjugate transpose

Eigenvalues and Eigenvectors

The eigenvalues and eigenvectors can be simply determined as well

e_val, e_vect = eig(a)

e_val will be an array of eigenvalues: e_vect is a 2-d array with column e_vect[:,i] containing the unit eigenvector corresponding to eigenvalue e_val[i].

Excercises

Please save the programs that you write to solve the exercises. We may ask you to show us the code and run it to re-generate the results.


1. Reading and analyzing 2 column data [1 ½ marks]

The file "data.txt" contains a two column set of data with a header. Edit your straight line fitting code to:

  1. Use the correct x and y labels, based on the header in the file. There is no need to write code to do this -- you can just open the file with a text editor, read the header and enter the values in your code manually.
  2. Read in the data from the rest of the file
  3. Use this information to make a plot of the data (with labelled axes) and add a best fit line to it.
  4. Create an output text file which contains the following information:
    • The names of the input data columns
    • The length of the data
    • Maxima and minima of the x and y data
    • The calculated intercept and gradient of the best fit line along with their errors.
    • The resistance in the circuit

Use string formatting techniques to write out an appropriate number of decimal places in your answers


2. Plotting 2-d arrays [2 marks]

Recall that the formula for the height of a wave with amplitude A, wavenumber k and angular frequency ω is y=A sin(kx-ωt) for a wave travelling to the right and y=A sin(kx+ωt) for a wave travelling to the left.

  1. Using mgrid(), generate 2-dimensional arrays for x and t both between 0 and 2 in steps of 0.01.
  2. Generate, a 2-dimensional array y giving the height of a wave with unit amplitude while the wavenumber and angular frequency are equal 2π.
  3. Make a contour plot using contourf() of the wave amplitude against x (on x-axis) and t (on y-axis).
  4. Make a 1-d plot of amplitude vs x at t=0, 0.25 and 0.5. (Hint:, it’s OK to use, e.g. t==0.25 to select the appropriate entries)
  5. Generate a standing wave by adding a second wave to y with same amplitude, wavenumber and frequency, but moving right.
  6. Make a contour plot of the standing wave amplitude against x and t.
  7. Make a 1-d plot of amplitude vs x at t=0, 0.25 and 0.5.

3. Matrix Manipulation [1 ½ marks]

Matrices A and B are given by:

$$A = \begin{pmatrix} 1.5 & 2\\ 0.5 & 1 \end{pmatrix} B = \begin{pmatrix} 3 & -2\\ -7 & 5 \end{pmatrix} $$

Write a program that asks the user for 8 numbers on one line. Split them, and put them into matrices A and B. The program should then calculate the following, and print out the answer.

  1. A-1
  2. B-1
  3. AB
  4. (AB)-1
  5. B-1 A-1
  6. The eigenvalues and eigenvectors for A

Note: we have not learned how to print matrices using string formatting from week 5 (i.e. the "{:.2f}" notation), so don't worry about this. Also, make sure you define A and B as matrices, not arrays.

When the demonstrator comes to mark this, they will put in different data than the A and B as defined above. The matrices above are reasonable "test" data to ensure your program works.