In [1]:
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
from __future__ import print_function
Jake VanderPlas
In this session, we'll discuss how to make your programs as efficient as possible, mainly by taking advantage of vectorization in NumPy.
If you use Python for any amount of time, you'll quickly find that there are some things it is not so good at. In particular, performing repeated operations via loops is one of its weaknesses.
For example:
In [2]:
# A silly function implemented in Python
def func_python(N):
d = 0.0
for i in range(N):
d += (i % 3 - 1) * i
return d
In [3]:
# Use IPython timeit magic to time the execution
%timeit func_python(10000)
To compare to a compiled language, let's write the same function in fortran and use the f2py tool (included in NumPy) to compile it
In [4]:
%%file func_fortran.f
subroutine func_fort(n, d)
integer, intent(in) :: n
double precision, intent(out) :: d
integer :: i
d = 0
do i = 0, n - 1
d = d + (mod(i, 3) - 1) * i
end do
end subroutine func_fort
In [5]:
!f2py -c func_fortran.f -m func_fortran > /dev/null
In [6]:
from func_fortran import func_fort
%timeit func_fort(10000)
Fortran is about 100 times faster for this task!
Languages tend to compromise between convenience and performance.
High level languages (Python, R, IDL, Matlab, etc.) tend to use High-level objects, dynamic typing, and interpreted execution for convenience in writing code
Low-level languages (C, Fortran, etc.) tend to use static types and compiled execution for performance
We like high-level languages because our development time is generally more valuable than execution time. But sometimes speed can be an issue.
Use Numpy ufuncs to your advantage
Use Numpy aggregates to your advantage
Use Numpy broadcasting to your advantage
Use Numpy slicing, masking, and fancy indexing to your advantage
Use a tool like SWIG, cython or f2py to interface to compiled code.
Here we'll cover the first four strategies.
A ufunc in numpy is a Universal Function. This is a function which operates element-wise on an array. We've already seen examples of these in the various arithmetic operations:
In [7]:
x = np.random.random(4)
print(x)
print(x + 1) # add 1 to each element of x
In [8]:
print(x * 2) # multiply each element of x by 2
In [9]:
print(x * x) # multiply each element of x by itself
In [10]:
print(x[1:] - x[:-1])
These are binary ufuncs: they take two arguments.
There are also many unary ufuncs:
In [11]:
-x
Out[11]:
In [12]:
np.sin(x)
Out[12]:
In [13]:
x = np.random.random(100000)
In [14]:
%%timeit
# compute element-wise x + 1 via a ufunc
y = x + 1
In [15]:
%%timeit
# compute element-wise x + 1 via a loop
y = np.empty_like(x)
for i in range(len(x)):
y[i] = x[i] + 1
Numpy UFuncs are faster than Python functions involving loops, because the looping happens in compiled code. This is only possible when types are known beforehand, which is why numpy arrays must be typed.
np.sin, np.cos, etc.)scipy.special.j0, scipy.special.gammaln, etc.)np.minimum, np.maximum)
In [16]:
x = np.random.random(5)
print(x)
print(np.minimum(x, 0.5))
print(np.maximum(x, 0.5))
In [17]:
# contrast this behavior with that of min() and max()
print(np.min(x))
print(np.max(x))
In [18]:
x = np.linspace(0, 10, 1000)
plt.plot(x, np.sin(x));
In [19]:
from scipy.special import gammaln
plt.plot(x, gammaln(x));
In [20]:
x = np.arange(5)
y = np.arange(1, 6)
x + y
Out[20]:
All operators (like +) actuall call an underlying numpy function: in this case np.add:
In [21]:
np.add(x, y)
Out[21]:
These ufuncs have some interesting and useful properties:
In [22]:
np.add.accumulate(x)
Out[22]:
In [23]:
np.multiply.accumulate(x)
Out[23]:
In [24]:
np.multiply.accumulate(y)
Out[24]:
In [25]:
np.add.identity
Out[25]:
In [26]:
np.multiply.identity
Out[26]:
In [27]:
np.add.outer(x, y)
Out[27]:
In [28]:
# make a times-table
x = np.arange(1, 13)
np.multiply.outer(x, x)
Out[28]:
Each of the following functions take an array as input, and return an array as output. They are implemented using loops, which is not very efficient.
For each function, implement a fast version which uses ufuncs to calculate the result more efficiently. Double-check that you get the same result for several different arrays.
use the %timeit magic to time the execution of the two implementations for a large array (say, 1000 elements).
In [29]:
# 1. computing the element-wise sine + cosine
from math import sin, cos
def slow_sincos(x):
"""x is a 1-dimensional array"""
y = np.zeros_like(x)
for i in range(len(x)):
y[i] = sin(x[i]) + cos(x[i])
return y
x = np.random.random(5)
slow_sincos(x)
Out[29]:
In [30]:
# write a fast_sincos function
In [31]:
# 2. computing the difference between adjacent squares
def slow_sqdiff(x):
"""x is a 1-dimensional array"""
y = np.zeros(len(x) - 1)
for i in range(len(y)):
y[i] = x[i + 1] ** 2 - x[i] ** 2
return y
x = np.random.random(5)
slow_sqdiff(x)
Out[31]:
In [32]:
# write a fast_sqdiff function
In [33]:
# 3. computing the outer-product of each consecutive pair
def slow_pairprod(x):
"""x is a 1-dimensional array"""
if len(x) % 2 != 0:
raise ValueError("length of x must be even")
N = len(x) / 2
y = np.zeros((N, N))
for i in range(N):
for j in range(N):
y[i, j] = x[2 * i] * x[2 * j + 1]
return y
x = np.arange(1, 9)
slow_pairprod(x)
Out[33]:
In [34]:
# write a fast_pairprod function
Aggregates are functions over arrays which return smaller arrays.
Numpy has several built-in
In [35]:
# 10 x 10 array drawn from a standard normal
x = np.random.randn(10, 10)
In [36]:
x.mean()
Out[36]:
In [37]:
np.mean(x)
Out[37]:
In [38]:
x.std()
Out[38]:
In [39]:
x.var()
Out[39]:
In [40]:
x.sum()
Out[40]:
In [41]:
x.prod()
Out[41]:
In [42]:
np.median(x)
Out[42]:
In [43]:
np.percentile(x, 50)
Out[43]:
In [44]:
np.percentile(x, (25, 75))
Out[44]:
In [45]:
x = np.random.rand(3, 5)
x
Out[45]:
In [46]:
x.sum(0) # sum along rows
Out[46]:
In [47]:
x.sum(1) # sum along columns
Out[47]:
In [48]:
np.median(x, 1)
Out[48]:
In [49]:
np.mean(x, 1)
Out[49]:
Any binary ufunc (a ufunc taking two arguments) can be turned into an aggregate using the reduce() method:
In [50]:
np.sum(x, 1)
Out[50]:
In [51]:
np.add.reduce(x, 1)
Out[51]:
In [52]:
np.prod(x, 1)
Out[52]:
In [53]:
np.multiply.reduce(x, 1)
Out[53]:
Some of these reduce methods are a bit silly:
In [54]:
np.divide.reduce(x, 1)
Out[54]:
A caution: for reduce methods, the default axis is 0, while for sum the default axis is None:
In [55]:
np.add.reduce(x)
Out[55]:
In [56]:
np.sum(x)
Out[56]:
Python has a min, max, and sum aggregate built-in. These are much more general than the versions in NumPy:
In [57]:
x = np.random.random(10000)
%timeit np.sum(x)
%timeit sum(x)
This is slow for the same reason that loops are slow: built-in min, sum, and max functions basically have to type-check every element.
Make sure to use Numpy's np.sum, np.min, and np.max.
Take the following functions, and convert them into an efficient form using aggregates. Each function expects a 1-dimensional array as input. Double-check that your function returns the same result as the original
In [58]:
def slow_cubesum(x):
"""x is a 1D array"""
result = 0
for i in range(len(x)):
result += x[i] ** 3
return result
x = np.random.random(100)
slow_cubesum(x)
Out[58]:
In [59]:
# implement fast_cubesum
In [60]:
def slow_rms(x):
"""x is a 1D array"""
m = np.mean(x)
rms = 0
for i in range(len(x)):
rms += (x[i] - m) ** 2
rms /= len(x)
return np.sqrt(rms)
x = np.random.random(100)
slow_rms(x)
Out[60]:
In [61]:
# implement fast_rms
Now we return to our silly function from the beginning of this section. Can you implement a fast version using ufuncs and aggregates?
In [62]:
def slow_sillyfunc(N):
"""N is an integer"""
d = 0.0
for i in range(N):
d += (i % 3 - 1) * i
return d
slow_sillyfunc(100)
Out[62]:
In [63]:
# Implement fast_sillyfunc using ufuncs & aggragates
Broadcasting is a way to extend the usefulness of ufuncs. For example, it allows you to subtract a vector from a matrix without writing a loop:
In [64]:
X = np.random.rand(3, 4)
mean = X.mean(0)
print(X)
print(mean)
In [65]:
print(X - mean)
We've subtracted the mean vector from each row of the matrix with a single operation. This is an example of broadcasting.
Broadcasting rules:
If the two arrays differ in their number of dimensions, the shape of the array with fewer dimensions is padded with ones on its leading (left) side.
If the shape of the two arrays does not match in any dimension, the array with shape equal to 1 in that dimension is stretched to match the other shape.
If in any dimension the sizes disagree and neither is equal to 1, an error is raised.
In [66]:
x = np.arange(10)
x ** 2
Out[66]:
In [67]:
Y = x * x[:, np.newaxis]
Y
Out[67]:
In [68]:
Y + 10 * x
Out[68]:
In [69]:
Y + 10 * x[:, np.newaxis]
Out[69]:
In [70]:
Y = np.random.random((2, 3, 4))
x = 10 * np.arange(3)
Y + x[:, np.newaxis]
Out[70]:
Now, assume you have $N$ points in $D$ dimensions, represented by an array of shape [N, D], try the following:
np.mean aggregate (that is, find the D-dimensional point which is the mean of the rest of the points)np.add ufunc.np.std aggregate.np.add ufunc.M, the centered and normalized version of the X array: $$ M_{ij} = (X_{ij} - \mu_j) / \sigma_j $$ This is one version of whitening the array.
In [71]:
X = np.random.random((1000, 5)) # 1000 points in 5 dimensions
In [72]:
# 1. Compute the mean of the 1000 points in X
In [73]:
# 2. Compute the mean using np.add
In [74]:
# 3. Compute the standard deviation across the 1000 points
In [75]:
# 4. Compute the standard deviation using np.add only
In [76]:
# 5. Compute the whitened version of the array
The last strategy we will cover is fancy indexing and masking.
For example, imagine you have an array of data where negative values indicate some kind of error.
In [77]:
x = np.array([1, 2, 3, -999, 2, 4, -999])
How might you clean this array, setting all negative values to, say, zero?
In [78]:
for i in range(len(x)):
if x[i] < 0:
x[i] = 0
x
Out[78]:
A faster way is to construct a boolean mask:
In [79]:
x = np.array([1, 2, 3, -999, 2, 4, -999])
mask = (x < 0)
mask
Out[79]:
And the mask can be used directly to set the value you desire:
In [80]:
x[mask] = 0
x
Out[80]:
Often you'll see this done in a single step:
In [81]:
x = np.array([1, 2, 3, -999, 2, 4, -999])
x[x < 0] = 0
x
Out[81]:
In [82]:
x = np.random.random(5)
x
Out[82]:
In [83]:
x[x > 0.5] = np.nan
x
Out[83]:
In [84]:
x[np.isnan(x)] = np.inf
x
Out[84]:
In [85]:
# warning: don't replace np.isnan(x) with (x == np.nan):
np.nan == np.nan
Out[85]:
In [86]:
x[np.isinf(x)] = 0
x
Out[86]:
In [87]:
x = np.array([1, 0, -1, -np.inf, np.inf, np.nan])
print("input ", x)
print("x < 0 ", x[x < 0])
print("x > 0 ", x[x > 0])
print("isinf ", x[np.isinf(x)])
print("isnan ", x[np.isnan(x)])
print("isposinf", x[np.isposinf(x)])
print("isneginf", x[np.isneginf(x)])
Boolean operators can be used on masks: order of operations can be a gotcha – be sure to use parentheses!
In [88]:
x = np.arange(16).reshape((4, 4))
x
Out[88]:
In [89]:
x[x < 5]
Out[89]:
In [90]:
x[~(x < 5)]
Out[90]:
In [91]:
x[(x < 10) & (x % 2 == 0)]
Out[91]:
In [92]:
x[(x > 3) & (x < 8)]
Out[92]:
Sum over a mask to find the number of True elements:
In [93]:
x = np.random.random(100)
print("array is length", len(x), "and has")
print((x > 0.5).sum(), "elements are greater than 0.5")
In [94]:
# clip is a useful function:
x = np.clip(x, 0.3, 0.6)
np.sum(x < 0.3), np.sum(x > 0.6)
Out[94]:
In [95]:
# works for 2D arrays as well
X = np.random.random((10, 10))
print((X < 0.1).sum())
In [96]:
x = np.random.random((3, 3))
x
Out[96]:
In [97]:
np.where(x < 0.5)
Out[97]:
In [98]:
x[x < 0.5]
Out[98]:
In [99]:
x[np.where(x < 0.5)]
Out[99]:
When you index with the result of a where function, you are using what is called fancy indexing: indexing with tuples
In [100]:
X = np.arange(16).reshape((4, 4))
X
Out[100]:
In [101]:
X[(0, 1), (1, 0)]
Out[101]:
In [102]:
X[range(4), range(4)]
Out[102]:
In [103]:
X.diagonal()
Out[103]:
In [104]:
# Note: diagonal() is read-only (though this will change in numpy 1.9)
X.diagonal() = 100
In [105]:
# Instead, try this:
X[range(4), range(4)] = 100
X
Out[105]:
In [106]:
X = np.arange(24).reshape((6, 4))
X
Out[106]:
In [107]:
i = np.arange(6)
np.random.shuffle(i)
i
Out[107]:
In [108]:
X[i] # X[i, :] is identical
Out[108]:
Fancy indexing also works for multi-dimensional index arrays
In [109]:
i2 = i.reshape(3, 2)
X[i2]
Out[109]:
In [110]:
X[i2].shape
Out[110]:
It's all about moving loops into compiled code:
Use Numpy ufuncs to your advantage (eliminate loops!)
Use Numpy aggregates to your advantage (eliminate loops!)
Use Numpy broadcasting to your advantage (eliminate loops!)
Use Numpy slicing and masking to your advantage (eliminate loops!)
Use a tool like SWIG, cython or f2py to interface to compiled code.