In [1]:
import numpy as np
In [2]:
np.random.seed(42) # Setting the random seed
In [3]:
# a vector: the argument to the array function is a Python list
v = np.random.rand(10)
v
Out[3]:
In [4]:
# a matrix: the argument to the array function is a nested Python list
M = np.random.rand(10, 2)
M
Out[4]:
We can index elements in an array using the square bracket and indices:
In [5]:
# v is a vector, and has only one dimension, taking one index
v[0]
Out[5]:
In [6]:
# M is a matrix, or a 2 dimensional array, taking two indices
M[1,1]
Out[6]:
If we omit an index of a multidimensional array it returns the whole row (or, in general, a N-1 dimensional array)
In [7]:
M[1]
Out[7]:
The same thing can be achieved with using : instead of an index:
In [8]:
M[1,:] # row 1
Out[8]:
In [9]:
M[:,1] # column 1
Out[9]:
We can assign new values to elements in an array using indexing:
In [10]:
M[0,0] = 1
In [11]:
M
Out[11]:
In [12]:
# also works for rows and columns
M[1,:] = 0
M[:,1] = -1
In [13]:
M
Out[13]:
Index slicing is the technical name for the syntax M[lower:upper:step] to extract part of an array:
In [14]:
a = np.array([1,2,3,4,5])
a
Out[14]:
In [15]:
a[1:3]
Out[15]:
Array slices are mutable: if they are assigned a new value the original array from which the slice was extracted is modified:
In [16]:
a[1:3] = [-2,-3]
a
Out[16]:
M[lower:upper:step]:
In [17]:
a[::] # lower, upper, step all take the default values
Out[17]:
In [18]:
a[::2] # step is 2, lower and upper defaults to the beginning and end of the array
Out[18]:
In [19]:
a[:3] # first three elements
Out[19]:
In [20]:
a[3:] # elements from index 3
Out[20]:
In [21]:
a = np.array([1,2,3,4,5])
In [22]:
a[-1] # the last element in the array
Out[22]:
In [23]:
a[-3:] # the last three elements
Out[23]:
In [24]:
A = np.array([[n+m*10 for n in range(5)]
for m in range(5)])
A
Out[24]:
In [25]:
# a block from the original array
A[1:4, 1:4]
Out[25]:
In [26]:
# strides
A[::2, ::2]
Out[26]:
Numpy arrays support two different way of storing data into memory, namely
The storage strategy is controlled by the parameter order of np.array
In [27]:
import numpy as np
FC = np.array([[1, 2, 3], [4, 5, 6],
[7, 8, 9], [10, 11, 12]], order='F')
In [28]:
CC = np.array([[1, 2, 3], [4, 5, 6],
[7, 8, 9], [10, 11, 12]], order='C')
In [29]:
FC[0, 1]
Out[29]:
In [30]:
CC[0, 1]
Out[30]:
In [31]:
FC.shape
Out[31]:
In [32]:
CC.shape
Out[32]:
Fancy indexing is the name for when an array or list is used in-place of an index:
In [33]:
row_indices = [1, 2, 3]
A[row_indices]
Out[33]:
In [34]:
col_indices = [1, 2, -1] # remember, index -1 means the last element
A[row_indices, col_indices]
Out[34]:
We can also index masks:
bool, then an element is selected (True) or not (False) depending on the value of the index mask at the position each element:
In [35]:
b = np.array([n for n in range(5)])
b
Out[35]:
In [36]:
row_mask = np.array([True, False, True, False, False])
b[row_mask]
Out[36]:
In [37]:
# same thing
row_mask = np.array([1,0,1,0,0], dtype=bool)
b[row_mask]
Out[37]:
This feature is very useful to conditionally select elements from an array, using for example comparison operators:
In [38]:
x = np.arange(0, 10, 0.5)
x
Out[38]:
In [39]:
mask = (5 < x)
mask
Out[39]:
In [40]:
x[mask]
Out[40]:
Alternatively, we can use the condition (mask) array directly within brackets to index the array
In [41]:
x[(5 < x)]
Out[41]:
Index slicing is the technical name for the syntax M[lower:upper:step] to extract part of an array
In [ ]:
In [ ]:
In [ ]:
Use fancy indexing to get all the elements of the previous matrix that are equals to 1
In [ ]:
In [ ]:
In [ ]:
Generate an array of numbers from 0 to 20 with step 0.5.
Extract all the values greater than a randomly generated number in the same range.
In [ ]:
Using function repeat, tile, vstack, hstack, and concatenate we can create larger vectors and matrices from smaller ones:
In [42]:
a = np.array([[1, 2], [3, 4]])
In [43]:
# repeat each element 3 times
np.repeat(a, 3)
Out[43]:
In [44]:
# tile the matrix 3 times
np.tile(a, 3)
Out[44]:
In [45]:
b = np.array([[5, 6]])
In [46]:
np.concatenate((a, b), axis=0)
Out[46]:
In [47]:
np.concatenate((a, b.T), axis=1)
Out[47]:
In [48]:
np.vstack((a,b))
Out[48]:
In [49]:
np.hstack((a,b.T))
Out[49]:
Generally, we want to avoid iterating over the elements of arrays whenever we can (at all costs). The reason is that in a interpreted language like Python (or MATLAB), iterations are really slow compared to vectorized operations.
However, sometimes iterations are unavoidable. For such cases, the Python for loop is the most convenient way to iterate over an array:
In [57]:
v = np.array([1,2,3,4])
for element in v:
print(element)
In [58]:
M = np.array([[1,2], [3,4]])
for row in M:
print("row", row)
for element in row:
print(element)
enumerate function to obtain both the element and its index in the for loop:
In [59]:
for row_idx, row in enumerate(M):
print("row_idx", row_idx, "row", row)
for col_idx, element in enumerate(row):
print("col_idx", col_idx, "element", element)
# update the matrix M: square each element
M[row_idx, col_idx] = element ** 2
In [60]:
# each element in M is now squared
M
Out[60]:
As mentioned several times by now, to get good performance we should try to avoid looping over elements in our vectors and matrices, and instead use vectorized algorithms. The first step in converting a scalar algorithm to a vectorized algorithm is to make sure that the functions we write work with vector inputs.
In [61]:
def Theta(x):
"""
Scalar implemenation of the Heaviside step function.
"""
if x >= 0:
return 1
else:
return 0
In [62]:
Theta(array([-3,-2,-1,0,1,2,3]))
OK, that didn't work because we didn't write the Theta function so that it can handle with vector input...
To get a vectorized version of Theta we can use the Numpy function np.vectorize. In many cases it can automatically vectorize a function:
In [63]:
Theta_vec = np.vectorize(Theta)
In [65]:
Theta_vec(np.array([-3,-2,-1,0,1,2,3]))
Out[65]:
Universal functions (Ufuncs) work on arrays, element-by-element, or on scalars.
Ufuncs accept a set of scalars as input, and produce a set of scalars as output.
In [66]:
Theta_ufunc = np.frompyfunc(Theta, 1, 1)
print("Result: ", Theta_ufunc(np.arange(4)))
Vectorizing code is the key to writing efficient numerical calculation with Python/Numpy.
That means that as much as possible of a program should be formulated in terms of matrix and vector operations, like matrix-matrix multiplication.
However, sometimes it may happen that we need to perform operations between ndarray objects of different size and shape so we might be tempted of performing this operation explicitly (i.e. explicit for loop).
Although we will get correct results (assumed we iterated arrays in the right way), this will produce very inefficient code.
For this reason, NumPy has a very important concept called Broadcasting.
Broadcasting is Numpy's terminology for performing mathematical operations between arrays with different shapes (assuming scalars being array with 0-dim).
In [72]:
D = np.array([ [0.3, 2.5, 3.5],
[2.9, 27.5, 0],
[0.4, 1.3, 23.9],
[14.4, 6, 2.3]])
And we also have another vector adj of values that contains some adjusting factors that we might want to apply to each sample (row) of data in D
In [73]:
adj = np.array([9, 4, 4])
In [79]:
%%timeit
# Create a new array filled with zeros, of the same shape as macros.
result = np.zeros_like(D)
# Now multiply each row of macros by cal_per_macro. In Numpy, `*` is
# element-wise multiplication between two arrays.
for i in range(D.shape[0]):
result[i, :] = D[i, :] * adj
result
This is a reasonable approach when coding in a low-level programming language: allocate the output, loop over input performing some operation, write result into output.
In Numpy, however, this is fairly bad for performance because the looping is done in (slow) Python code instead of internally by Numpy in (fast) C code
(we've proven this already when we compared ndarray vs list)
In [80]:
%%timeit
adj_stretch = np.tile(adj, (D.shape[0], 1))
D * adj_stretch
In [81]:
adj_stretch
Out[81]:
Nice, it's shorter too, and slightly faster! To appreciate even more performance gain, of our np.tile solution, we could try increasing the size of D to a bigger structure:
In [84]:
D_large = np.random.rand(10**6, 10)
adj_large = np.random.rand(10)
In [85]:
D_large.shape, adj_large.shape
Out[85]:
In [86]:
%%timeit
# Create a new array filled with zeros, of the same shape as macros.
result_large = np.zeros_like(D_large)
# Now multiply each row of macros by cal_per_macro. In Numpy, `*` is
# element-wise multiplication between two arrays.
for i in range(D_large.shape[0]):
result_large[i, :] = D_large[i, :] * adj_large
In [87]:
%%timeit
adj_large_stretch = np.tile(adj_large, (D_large.shape[0], 1))
D_large * adj_large_stretch
The loop-in-Python method takes ~1.5 seconds, the stretching method takes ~48 milliseconds, a ~75x speedup.
And now, finally, comes the interesting part.
You see, the operation we just did - stretching one array so that its shape matches that of another and then applying some element-wise operation between them - is actually pretty common.
This often happens when we want to take a lower-dimensional array and use it to perform a computation along some axis of a higher-dimensional array.
In fact, when taken to the extreme this is exactly what happens when we perform an operation between an array and a scalar - the scalar is stretched across the whole array so that the element-wise operation gets the same scalar value for each element it computes.
Numpy generalizes this concept into broadcasting - a set of rules that permit element-wise computations between arrays of different shapes, as long as some constraints apply.
Incidentally, this lets Numpy achieve two separate goals - usefulness as well as more consistent and general semantics.
Binary operators like * are element-wise, but what happens when we apply them between arrays of different shapes? Should it work or should it be rejected? If it works, how should it work?
Broadcasting defines the semantics of these operations.
In [89]:
%%timeit
## Back to our example
D * adj # Broadcasting
Element-wise operations on arrays are only valid when the arrays' shapes are either equal or compatible.
The equal shapes case is trivial - this is the stretched array from the example above. What does compatible mean, though?
To determine if two shapes are compatible, Numpy compares their dimensions, starting with the trailing ones and working its way backwards.
For example, for the shape (4, 3, 2) the trailing dimension is 2, and working from 2 "backwards" produces: 2, then 3, then 4.
If two dimensions are equal, or if one of them equals 1, the comparison continues. Otherwise, you'll see a ValueError raised (saying something like "operands could not be broadcast together with shapes ...").
When one of the shapes runs out of dimensions (because it has less dimensions than the other shape), Numpy will use 1 in the comparison process until the other shape's dimensions run out as well.
Once Numpy determines that two shapes are compatible, the shape of the result is simply the maximum of the two shapes' sizes in each dimension.
More here: Broadcasting Arrays in NumPy
We can use the usual arithmetic operators to multiply, add, subtract, and divide arrays with scalar numbers.
In [90]:
v1 = np.arange(0, 5)
In [91]:
v1 * 2
Out[91]:
In [92]:
v1 + 2
Out[92]:
In [93]:
A = np.array([[n+m*10 for n in range(5)] for m in range(5)])
In [94]:
print('A * 2: ', '\n', A * 2)
print('A + 2: ', '\n', A + 2)
When we add, subtract, multiply and divide arrays with each other, the default behaviour is element-wise operations:
In [95]:
A * A # element-wise multiplication
Out[95]:
In [96]:
v1 * v1
Out[96]:
If we multiply arrays with compatible shapes, we get an element-wise multiplication of each row:
In [97]:
A.shape, v1.shape
Out[97]:
In [99]:
A * v1 #Broadcasting
Out[99]:
What about matrix mutiplication?
There are two ways.
We can either use the np.dot function, which applies a matrix-matrix, matrix-vector, or inner vector multiplication to its two arguments:
In [100]:
np.dot(A, A)
Out[100]:
In [101]:
np.dot(A, v1)
Out[101]:
In [102]:
np.dot(v1, v1)
Out[102]:
In [103]:
A @ v1
Out[103]: