.. For doctests >>> import numpy as np >>> np.random.seed(0) >>> # For doctest on headless environments >>> import matplotlib.pyplot as plt >>> plt.switch_backend("Agg") .. currentmodule:: numpy Numerical operations on arrays ============================== .. contents:: Section contents :local: :depth: 1 Elementwise operations ---------------------- Basic operations ................ With scalars: .. sourcecode:: pycon >>> a = np.array([1, 2, 3, 4]) >>> a + 1 array([2, 3, 4, 5]) >>> 2**a array([ 2, 4, 8, 16]) All arithmetic operates elementwise: .. sourcecode:: pycon >>> b = np.ones(4) + 1 >>> a - b array([-1., 0., 1., 2.]) >>> a * b array([2., 4., 6., 8.]) >>> j = np.arange(5) >>> 2**(j + 1) - j array([ 2, 3, 6, 13, 28]) These operations are of course much faster than if you did them in pure python: .. sourcecode:: pycon >>> a = np.arange(10000) >>> %timeit a + 1 # doctest: +SKIP 10000 loops, best of 3: 24.3 us per loop >>> l = range(10000) >>> %timeit [i+1 for i in l] # doctest: +SKIP 1000 loops, best of 3: 861 us per loop .. warning:: **Array multiplication is not matrix multiplication:** .. sourcecode:: pycon >>> c = np.ones((3, 3)) >>> c * c # NOT matrix multiplication! array([[1., 1., 1.], [1., 1., 1.], [1., 1., 1.]]) .. note:: **Matrix multiplication:** .. sourcecode:: pycon >>> c.dot(c) array([[3., 3., 3.], [3., 3., 3.], [3., 3., 3.]]) .. topic:: **Exercise: Elementwise operations** :class: green * Try simple arithmetic elementwise operations: add even elements with odd elements * Time them against their pure python counterparts using %timeit. * Generate: * [2**0, 2**1, 2**2, 2**3, 2**4] * a_j = 2^(3*j) - j Other operations ................ **Comparisons:** .. sourcecode:: pycon >>> a = np.array([1, 2, 3, 4]) >>> b = np.array([4, 2, 2, 4]) >>> a == b array([False, True, False, True]) >>> a > b array([False, False, True, False]) .. tip:: Array-wise comparisons: .. sourcecode:: pycon >>> a = np.array([1, 2, 3, 4]) >>> b = np.array([4, 2, 2, 4]) >>> c = np.array([1, 2, 3, 4]) >>> np.array_equal(a, b) False >>> np.array_equal(a, c) True **Logical operations:** .. sourcecode:: pycon >>> a = np.array([1, 1, 0, 0], dtype=bool) >>> b = np.array([1, 0, 1, 0], dtype=bool) >>> np.logical_or(a, b) array([ True, True, True, False]) >>> np.logical_and(a, b) array([ True, False, False, False]) **Transcendental functions:** .. sourcecode:: pycon >>> a = np.arange(5) >>> np.sin(a) array([ 0. , 0.84147098, 0.90929743, 0.14112001, -0.7568025 ]) >>> np.log(a) array([ -inf, 0. , 0.69314718, 1.09861229, 1.38629436]) >>> np.exp(a) array([ 1. , 2.71828183, 7.3890561 , 20.08553692, 54.59815003]) **Shape mismatches** .. sourcecode:: pycon >>> a = np.arange(4) >>> a + np.array([1, 2]) # doctest: +SKIP Traceback (most recent call last): File "", line 1, in ValueError: operands could not be broadcast together with shapes (4) (2) *Broadcasting?* We'll return to that :ref:later . **Transposition:** .. sourcecode:: pycon >>> a = np.triu(np.ones((3, 3)), 1) # see help(np.triu) >>> a array([[0., 1., 1.], [0., 0., 1.], [0., 0., 0.]]) >>> a.T array([[0., 0., 0.], [1., 0., 0.], [1., 1., 0.]]) .. note:: **The transposition is a view** The transpose returns a *view* of the original array:: >>> a = np.arange(9).reshape(3, 3) >>> a.T[0, 2] = 999 >>> a.T array([[ 0, 3, 999], [ 1, 4, 7], [ 2, 5, 8]]) >>> a array([[ 0, 1, 2], [ 3, 4, 5], [999, 7, 8]]) .. note:: **Linear algebra** The sub-module :mod:numpy.linalg implements basic linear algebra, such as solving linear systems, singular value decomposition, etc. However, it is not guaranteed to be compiled using efficient routines, and thus we recommend the use of :mod:scipy.linalg, as detailed in section :ref:scipy_linalg .. topic:: Exercise other operations :class: green * Look at the help for np.allclose. When might this be useful? * Look at the help for np.triu and np.tril. Basic reductions ---------------- Computing sums .............. .. sourcecode:: pycon >>> x = np.array([1, 2, 3, 4]) >>> np.sum(x) 10 >>> x.sum() 10 .. image:: images/reductions.png :align: right Sum by rows and by columns: .. sourcecode:: pycon >>> x = np.array([[1, 1], [2, 2]]) >>> x array([[1, 1], [2, 2]]) >>> x.sum(axis=0) # columns (first dimension) array([3, 3]) >>> x[:, 0].sum(), x[:, 1].sum() (3, 3) >>> x.sum(axis=1) # rows (second dimension) array([2, 4]) >>> x[0, :].sum(), x[1, :].sum() (2, 4) .. tip:: Same idea in higher dimensions: .. sourcecode:: pycon >>> x = np.random.rand(2, 2, 2) >>> x.sum(axis=2)[0, 1] # doctest: +ELLIPSIS 1.14764... >>> x[0, 1, :].sum() # doctest: +ELLIPSIS 1.14764... Other reductions ................ --- works the same way (and take axis=) **Extrema:** .. sourcecode:: pycon >>> x = np.array([1, 3, 2]) >>> x.min() 1 >>> x.max() 3 >>> x.argmin() # index of minimum 0 >>> x.argmax() # index of maximum 1 **Logical operations:** .. sourcecode:: pycon >>> np.all([True, True, False]) False >>> np.any([True, True, False]) True .. note:: Can be used for array comparisons: .. sourcecode:: pycon >>> a = np.zeros((100, 100)) >>> np.any(a != 0) False >>> np.all(a == a) True >>> a = np.array([1, 2, 3, 2]) >>> b = np.array([2, 2, 3, 2]) >>> c = np.array([6, 4, 4, 5]) >>> ((a <= b) & (b <= c)).all() True **Statistics:** .. sourcecode:: pycon >>> x = np.array([1, 2, 3, 1]) >>> y = np.array([[1, 2, 3], [5, 6, 1]]) >>> x.mean() 1.75 >>> np.median(x) 1.5 >>> np.median(y, axis=-1) # last axis array([2., 5.]) >>> x.std() # full population standard dev. 0.82915619758884995 ... and many more (best to learn as you go). .. topic:: **Exercise: Reductions** :class: green * Given there is a sum, what other function might you expect to see? * What is the difference between sum and cumsum? .. topic:: Worked Example: diffusion using a random walk algorithm .. image:: random_walk.png :align: center .. tip:: Let us consider a simple 1D random walk process: at each time step a walker jumps right or left with equal probability. We are interested in finding the typical distance from the origin of a random walker after t left or right jumps? We are going to simulate many "walkers" to find this law, and we are going to do so using array computing tricks: we are going to create a 2D array with the "stories" (each walker has a story) in one direction, and the time in the other: .. only:: latex .. image:: random_walk_schema.png :align: center .. only:: html .. image:: random_walk_schema.png :align: center :width: 100% .. sourcecode:: pycon >>> n_stories = 1000 # number of walkers >>> t_max = 200 # time during which we follow the walker We randomly choose all the steps 1 or -1 of the walk: .. sourcecode:: pycon >>> t = np.arange(t_max) >>> steps = 2 * np.random.randint(0, 1 + 1, (n_stories, t_max)) - 1 # +1 because the high value is exclusive >>> np.unique(steps) # Verification: all steps are 1 or -1 array([-1, 1]) We build the walks by summing steps along the time: .. sourcecode:: pycon >>> positions = np.cumsum(steps, axis=1) # axis = 1: dimension of time >>> sq_distance = positions**2 We get the mean in the axis of the stories: .. sourcecode:: pycon >>> mean_sq_distance = np.mean(sq_distance, axis=0) Plot the results: .. sourcecode:: pycon >>> plt.figure(figsize=(4, 3)) # doctest: +ELLIPSIS