Arrays are the workhorse for engineers and scientists.
Arrays are part of the numpy package. For a good overvies see here.
For visualization we use matplotlib as always.
In [272]:
import matplotlib.pyplot as plt
import numpy as np
"""You may have to uncomment one of the two following lines to make in line plots available.
But this is not always necessary, for instance not on the mac.
"""
#%matplotlib inline
#%matplotlib notebook # for intractive plots
Out[272]:
Arrays
Arrays are contiguous blocks of data in memory with all elements of the same type. That is mostly ints or floats, but can be anything. This type is called dtype. Arrays can have as many dimensions as you like and may be very large. THey can be treated as a single unit, which makes them very powerful.
There are many ways to create arrays. One is to type them.
A 1D array is looks a bit like a list but it is not. A 2D array looks like a list of lists. A 3D array looks like a list of lists of lists and so on.
In [273]:
# A one-D array
A = np.array([1, 7, 2, 12]) # specify a one-dimensional array
# A 2D array
B = np.array([[1, 3, 5], [3, 1, 6.2], [4.1, 6, 1], [4, 1, 3]])
# A 2D array
C = np.array([[ 1, 3, 5],
[ 3, 1, 6.2],
[ 4.1, 6, 1],
[ 4, 1, 3]])
#A 3D array
D = np.array([ [ [ 1, 2, 3],
[ 4, 5, 6],
[ 7, 8, 9] ],
[ [10, 11, 12],
[13, 14, 15],
[16, 17, 18] ],
[ [19, 20, 21],
[22, 23, 24],
[25, 26, 27] ] ])
print("A = \n", A, "\n")
print("B = \n", B, "\n")
print("C = \n", C, "\n")
print("D = \b", D, "\n")
In [237]:
print("The type of the data in array A is :", A.dtype)
print("The type of the data in array C is :", C.dtype)
print("The type of the data in array D is :", D.dtype)
print("Size of A: ", A.size)
print("size of C: ", C.size)
print("Size of D: ", D.size)
print("Shape of A: ", A.shape, " <-- Notice the comma !")
print("Shape of C: ", C.shape)
print("Shape of D: ", D.shape)
print("Nr of dimensions of A: ", A.ndim)
print("Nr of dimensions of C: ", C.ndim)
print("Nr of dimensions of D: ", D.ndim)
In [238]:
A = np.array([1, 2, 3, 4], dtype=float)
print("Array A :", A, " now has dtype=", A.dtype)
To fast generate arrays of ones or zeros we can use the function np.ones(tuple) or np.zeros(tuple).
In [239]:
A = np.zeros((4, 5)) # note (4, 5) is a tuple, parantheses needed.
B = np.ones((3, 6))
print("A =\n", A)
print("B =\n", B)
Or generate an array with a shape that is the same as an already defined other array:
In [240]:
C = np.ones_like(A)
print("C =\n", C)
Use the function np.arange(start, end, step) to generate a regular sequence.
In [241]:
print(np.arange(1, 6)) # start, end with default steps 1
print(np.arange(7)) # starts at 0 end ends at 6 (7 numbers in total)
print(np.arange(2, 20, 3))
A practical way to generate one-dimensional arrays is by using the functions linspace(start, end, number of points)
and logspace(log10(start), log10(end), number of points)
In [242]:
x = np.linspace(0, 6 * np.pi, 201)
y = np.sin(x)
z = y * np.exp(-x/5.0)
plt.figure()
plt.plot(x, y, label= '$sin(x)$')
plt.plot(x, z, label= '$sin(x) e^{-x/(2 \pi)}$')
plt.title('Some graphs')
plt.xlabel('x')
plt.ylabel('y')
plt.legend(loc='best', fontsize='small')
plt.show()
Individual elements of an array can be addressed by index. Indices in Python start at 0 and are placed in square brackets.
Some examples
In [243]:
x = np.arange(20, 30) # starts at 20 and stops before 30, 10 elements
print(x)
print(x[0:5]) # 0 through 4 is 5 elements
print(x[:5]) # same, zero assumed by Python
print(x[3:7]) # 3 through 6, 7-3 = 4 elements
print(x[2:9:2]) # starts at 2 stops before 9 with steps equal to 2
Backcounting:
In [244]:
print(x[-1:0:-1]) # starts at last ends before first with steps -1
print(x[-1:0:-2]) # starts at last end before first with steps -2
print(x[-1::-1]) # starts at last ends at first with setps of -1 !! (note the difference)
More-dimensional arrays. It's easy to generate multidimensional arrays from one-dimensional ones using the reshape function. For instance to show the element number of an array at which it is stored in memory:
In [245]:
A = np.arange(3 * 4 * 5).reshape((3, 4, 5))
print("3D array A: \n", A)
The easiest way to work with 3D spatial arrays is to consider them as [Layer Row Col] or [z, y, x]. When looked upon in this way, each bock is a layer, each row is along the y-axis on paper ane each row along the x-axis.
Of course, you could associate the first dimension with y the second with x and the thrid with z, as it s natural in Matlab, but that is completely confusing in Python. This is because the elements in a numpy array are stored with the last dimension changing fastest (Fortran order), while in Matlab it's the other way around (C order). This has consequences for the way arrays are printed. To match the print with the imagination, stick with the layer, row, col association in Python numpy arrays.
The reshape function often very useful. Notice that it requires a tuple with the dimensions and it requires that the number of elements in the array does not change by the reshape.
We can address any cell by using indices:
In [246]:
print(A[2, 2, 1])
This is is layer 3, row 3 column 2 (because indexing starts at 0).
In [247]:
print(A[1:3, 0:2, 1:3])
That is: of layer 1 and 2 we take rows 0, and 1 and of them columns 1 and 2.
The other way around, starting at the last index: Columns 1:3 (= 1 and 2) of rows 0:2 (=0 and 1) of layers 1:3 (= 1 and 2)
We can simply select layers by using only the first dimension:
In [248]:
print(A[1:3]) # selects block 1 and 2
And the rows of them by using only the first 2 dimensions:
In [249]:
print(A[1:3, 0:2])
And, finally
In [250]:
print(A[1:3, 0:2, 1:3])
You can assign one value to a range of an array by specifying a range of indices,
or you can assign an array to a range of another array, as long as the ranges have equal length. In the last example below, the first 5 values of x (specified as x[0:5]) are given the values [40,42,44,46,48].
In [251]:
print(A[1:3])
Indexing in sequence has a different meaning:
In [252]:
print(A[1:3][0:2][1:3])
This should be read as
In [253]:
(((A[1:3])[0:2])[1:3])
Out[253]:
A[1:3] picks blocks 1 and 2 of the 3D array A, which is also a 3D array. Of that, the next index [0:2] takes blocks 0 and 1, which is, in fact the entire array A[1:3] Of that array, the index [1:3] picks blocks 1 and 2.
Notice that if the end is larger than the size of the array in that dimensions, its value defaults to the last. In some languages this would be considered a bug.
In [254]:
((A[1:3])[1:10000])[0:1000]
Out[254]:
It's easy to revert an array:
In [255]:
A[-1::-1, -1::-1, -1::-1]
Out[255]:
Or with the well-known short hand:
In [256]:
A[::-1, :, ::-1] # only the layers and columns reverted, not the rows
Out[256]:
Arrays use one dimension if only one layer, row or column in selected:
In [257]:
print("A=\n", A)
print("\nA.shape = ", A.shape, " <--- 3 dimensions")
print("\nA.ndim = ", A.ndim)
print("\nA[1]=\n", A[1])
print("\nA[1].shape =", A[1].shape, " <--- 2 dimensions")
print("\nA[1].ndim =", A[1].ndim)
print("\nA[1][1] =\n", A[1][1])
print("\nA[1][1].shape =", A[1][1].shape, " <--- 1 dimension")
print("\nA[1][1].ndim =", A[1][1].ndim)
print("\nA[1][1][1] =", A[1][1][1])
print("\nA[1][1][1].shape =", A[1][1][1].shape, "<--- 0 dimensions")
print("\nA[1][1][1].ndim =", A[1][1][1].ndim)
Broadcasting
Arrays of different size can be broadcasted if the broadcasting rul applies.
The broadcasting rule
Two arrays are compatible for broadcasting if that is, starting from the end, the axis lengths match or if either of the lengths is 1. Broadcasting is then performed over the missing and / or length 1 dimensions.
This sounds very complicated, but one may also say. Both arrays must have the same number of dimensions and the lengths along the dimimsions in the two arrray must either match or be 1.
For example:
We want to generate a 3D array with conductivities, one for every layer as defined by a vector. The array in the xy plane has 5 by 4 elements
In [258]:
D = np.ones((5, 4)) # size of array in xy plane
k = np.array([4.2, 11.7, 13.1, 25.2, 33.]) # layer conductivities
To generate the desired array K, we have to turn the vector k in to a 3D vector
In [259]:
print("Original k.shape :", k.shape, " <-- 1D")
k = k.reshape((len(k), 1, 1))
print("New k.shape :", k.shape, " <-- 3D")
Still we cannot multiply by D because D is not threedimensional. So turn D into a 3D array with only one layer
In [260]:
print("Origianal D.shape :", D.shape) # <-- 2D
D = D.reshape((1, *D.shape))
print("New D.shape :", D.shape) # <-- 3D
Now generate the conductivity array K
In [261]:
K = D * k
print("K =\n", K)
A few more examples to exercise
First, only one dimension has length 1.
In [262]:
n, m, k = A.shape
B = np.ones((1, m, k))
A * B
C = np.ones((n, 1, k))
A * C
D = np.ones((n, m, 1))
A * D
A * B * C * D
Out[262]:
Next two dimensions have lengt 1.
In [263]:
n, m, k = A.shape
B = 3 * np.ones((n, 1, 1))
A * B
C = 3 * np.ones((1, m, 1))
A * C
D = 3 * np.ones((1, 1, k))
A * D
A * B * C * D
Out[263]:
B, C an D are vectors pointing in the z, y an x directions respectively:
In [264]:
n, m, k = A.shape
B = A[:, 0:1, 0:1]
print("B=\n", B)
C = A[0:1, :, 0:1]
print("C=\n", C)
D = A[0:1, 0:1, :]
print("D=\n", D)
print("B * C * D = \n", B * C * D)
In [265]:
B = np.random.randint(0, 17, size=(7, 3))
print("B = \n", B)
print("\nB*B=\n", B * B)
print("\nnp.dot(B.T, B)=\n", np.dot(B.T, B))
print("\n B @ B = \n", B.T @ B) # is the same as np.dot(B.T, B)
For true linear algebra, one may choose to use the np.linalg package which provides matrices and linear algebra functions.
With the array converted to a matrix the multiplication is treated as a matrix multiplication.
In [266]:
print("Using linear algebra, matrix:")
B = np.matrix(B)
print("\nB.T * B = \n ", B.T * B)
print("\nInverse of B = B.I = \n",B.I)
Using conditional indexing
A condition on an array like y<20 yields a boolean array of the same dimension as y with the value True wherever the condition is met and False otherwise. We can use conditions, i.e. boolean arrays, on for logical indexing of arrays. Only the values where the boolean array is True will be used.
so if y has 100 values and y<3 is true for only 20 values, then y[y<3] also has only 20 values because only the y's are picked from the array y where the condition is True. Logical indexing is very often the most effective, intuitive and elegant way of picking parts from arrays, without using any numerical indices.
In [271]:
x = np.linspace(0, 4 * np.pi, 201)
y = np.sin(x)
z = y * np.exp(-x/5.)
w = (x - 2.3) * (x - 5.1 ) * (x - 10.3) * (x - 12)
# some logical indexing
w[w > z] = z[w > z]
w[w < y] = y[w < y]
plt.plot(x, y, label="y")
plt.plot(x, z, label="z")
plt.plot(x, w, label="w", lw=2)
plt.legend()
plt.show()
# or, to show how intuitive this indexing is:
high = w > z # w > z yields a boolean array
low = w < y
w[high] = z[high]
w[low ] = y[low]
plt.plot(x, w, label="w", lw=2)
plt.show()
In [269]:
AND = np.logical_and
x = np.linspace(-500, 500, 200)
y = np.linspace(-500, 500, 200)
xBox = [-300, 50]
yBox = [-200, 300]
R0 = 300.
x0, y0 = 50, 120
y = y.reshape((len(y), 1))
x = x.reshape((1, len(x)))
X = np.ones_like(y) * x
Y = y * np.ones_like(x)
R = np.sqrt((X - x0)**2 + (Y - y0)**2)
inBox = AND(AND(X > xBox[0], X < xBox[1]), AND(Y > yBox[0], Y < yBox[1]))
inIsland = R < R0
V = np.zeros_like(X)
V[inBox] = 1
V[inIsland] = 2
plt.matshow(V)
plt.show()