CHAPTER NumPy: creating and manipulating numerical data Authors: Emmanuelle Gouillart, Didrik Pinte, Gaél Varoquaue, and Pauli Virtanen This chapter gives an overview of NumPy, the core yerformant numerical computing with Python, 4.1 The NumPy array object Section contents hat are NumPy and NumPy arrays? 4.1.1 What are NumPy and NumPy arrays? NumPy arrays Python objects 50Scipy lecture notes, Edi + high-level number objects: integers, floating point + containers: lists (costless insertion and append), dictionaries (fast lookup) NumPy provides + extension package to Python for multi-dimensional arrays + closer to hardware (efficiency) + designed for scientific computation (convenience) + Also known as array oriented computing >>> import mumpy as =p >>> a= mp.array([0, 1, 2, 31) array({0, 1, 2, 31) ‘Tip: For example, An array containing + values of an experiment/sinaulation at discrete time steps + signal recorded by a measurement device, e.g, sound wave + pixels of an image, grey-level or colour + 3:D data measured at different X-Y-Z positions, e.g. MRI scan Why it is useful: Memory-efficient container that provides fast numerical operations. Im [1]: b= range(1000) Im (2): feimest [1442 for a an LI 1000 loeps, best of 3: 403 us per loop Im [3]: a = mp.arange(1000) Im [4]: feimest as2 100000 loops, best of 3: NumPy Reference documentation + On the web: bttps://numpy:org/doc/ + Interactive help: In (5): np.arrayf] String Forn:cbuilt-in function array> Docstring: array (object, dtyp fone, copy=True, ode! ne, subok-False, nénis + Looking for something: >>> mp leokfor('create array") Search results for ‘ere array! umpy. array (conties om He Dae) 4.1, The NumPy array object 51Scipy lecture notes, Edition 2022.1 {continued from previous page) ore rumpy.menaap an array ‘a nenory-map to an array stored in a woinarys file on disk Ta [6]: mp-com-B] np.conj np conjugate np.convolve Import conventions ‘The recommended convention to import mumpy is: >>> import munpy ae =p 4.1.2 Creating arrays Manual construction of arrays +L de> a= mp.array((0, 1, 2, 31) array([0, 1, 2, 31) >>> a.ndim. 1 >>> a.shape a >>> lena) a + 2D, 3D,... >>> b= mp.array(([0, 1, 23, [3, 4, 1) #225 array >>> array({(0, 1, 2] (3, 4, sl) >> binaia, 2 D> shape 2, 3) >>> len(b) returns the stze of the first dimension 2 pope = mp-array((0[i], (211, (131, (410) array(C((1), a. ta), tan >>> ¢.shape 4.1, The NumPy array object 52Scipy lecture notes, Edition 2022.1 + Greate a simple two dimensional array. First, redo the examples from above. And then create your own: how about odd numbers counting backwards on the first row, and even numbers on the second? + Use the funetions 1en(), numpy.shape() on these arrays. How do they relate to each other? And to the ndim attribute of the arrays? Functions for creating arrays ‘Tip: In practice, we rarely enter items one by oue. + Evenly spaced: array([0, 1, 2, 3, 4, 5, 6, 7, 8, 91) >>> b= mp.arange(1, 9, 2) # Start, end (exclusive), step >>> array((1, 3, 5, 71) ap-arange(10) #0... m1 (1) + of by number of points: 55> © = np. lanepace(0, array([0. , 0.2, 0.4, >>> d= np Linspace(0, >> array(l0. , 0.2, 0.4, + Common arrays: >>> a= mp.ones((3, 3)) # reminder: (3, 3) ts o tuple array({Li., 4., 11, Gy 4, 42, fy 4, 1p de> b= np zeroe((2, 2)) >>> b array({[0., 0-1, [o., 0.11) >>> ¢ = up_eye(3) array({[t., 0., 0.1, fo, 1, 0.1, fo, 0, 1.10) de> d= np-diaginp.array((1, 2, 3, 41)) aod + mp.xandom: random numbers (Mersenne Twister PRNG) >>> a= np. random. rand(4) ¥ uniform én (0, 1) array({ 0.95799151, 0.14222247, 0.08777354, 0.51887988)) >> b np-randon.randn(s) # Gaussian (Gontinaes om we BRR) 4.1, The NumPy array object 53Scipy lecture notes, Edition 2022.1 {continued from previous page) >> be array({ 0.37644698, 0. 11426368, -0,47616538, 1,79664113)) >>> mp.randon.seed(1234) # Setting the random seed Exercise: Creating arrays using functions + Experiment with arange, linspace, ones, zeros, eye and diag. + Create different kinds of arrays with random numbers, + Try setting the seed before ereating an array with random values. + Look at the function np.empty. What does it do? When might this be useful? 4.1.3 Basic data types You may have noticed that, in some instances, array elements are displayed with a trailing dot (e.g. 2 vs 2). This is due to a difference in the data-type used: >o> a= mp-erray([1, 2, 3]) >> a.dtype atype(*ines4*) de> b= mp.array(L1., 2, 3.1) >>> b.dtype type ('sloac64") Different data-types allow us to store data more compactly in memory, but most of the time we uply work with floating point numbers. Note that, in the example above, NamPy auto-detects the data-type from the input, You can explicitly specify which data-type you want: de> e = mp.array([1, 2, 3], atyperfloat) >>> ©. dtype type ('sloac64") The default data type is floating point 55> a= mp.ones((3, 3)) >>> a.dtype atype('sleac64") ‘There are also other types: Complex >>> d= mp.array([i+2j, 3+4), 5+542]]) >>> a.atype type ( conplext28') Bool >>> @ = np.array([True, False >>> e.atype type ("bool") False, True)) 4.1, The NumPy array object 54Scipy lecture notes, Edition 2022.1 Strings >>> f= mp.array(['Bonjour', ‘Helle’, ‘Halle']) >>> fdtype —# <--- strings containing maz. 7 letters avype('S7") Much more + anta2 + anted + wines + wines 4.1.4 Basic visualization Now that we have our first data arrays, we are going to visualize them. Start by launching IPython: $ apython # or épython3 depending on your install Or the notebook: $ jupyter notebook Once Python has started, enable interactive plots: 55> Ymatploeiib Or, from the notebook, enable plots in the notebook: >>> lmatplotiib inline ‘The inline is important for the notebook, so that plots are displayed in the notebook and not in a new window Matplotlib is a 2D plotting package, We can import its functions as below: >>> import matplotlib.pyplot as plt # the tidy woy And then use (note that you have to use show explicitly if you have not enabled interactive plots with ‘imatplotib) >>> plt.plec(x, y) # Line plot >>> plt-shew() # <= shows the plot (not needed with interactive plots) Or, if you have enabled interactive plots with Yjmatplot1ib: >>> ple plos(x, y) # Line plot + ID plotting: >>> x = mp. Linspace(0, 3, 20) >>> y = np-Lingpace(0, 9, 20) >>> plt-plestx, y) # line plot [cnatplotiib-Lines.Line2D object at ...>) de> plt.plot(x, y, '0') # dot plot [enatplotisb. lines Line2D object at ...>] 4.1. The NumPy array object 55Scipy lecture notes, Editi + 2D arrays (such as images) >>> image = p.randem.rand(30, 30) de> pt. imshow(inage, cnap=plt.cm.hot) cnatplotlib. image. xestnage object at >>> ple. colerbar() >> a= mp-arange(10) array((9, 1, 2, 3, 4, 5, 6, 7%, 8 91) >>> alo], al2], al-1] (0, 2, 9) 4.1, The NumPy array object 56Scipy lecture notes, Edition 2022.1 Warning: Indices begin at 0, like other Python sequences (and C/C++). In contrast, in Fortran, or Matlab, indices begin at 1. ‘The usual python idiom for reversing a sequence is supported: be alii] array([9, 8, 7, 6, 5, 4, 3, 2, 1, OD For multidimensional arrays, indices are tuples of integers’ >>> a= np. diag(mp. arange(3)) array({{0, 0, 0}, 0, 1, 0), [o, 0, 211) >>> alt, 12 1 >>> al2, 1] = 10 # third Line, second colum array({ 0, 0, ol, Lo, 1, ol, Lo, 10, 21) >> alt] array([0, 1, 0) Note: + In 2D, the first dimension corresponds to rows, the second to columns. + for multidimensional a, a[0] is interpreted by taking all clements in the unspecified dimensions. Slicing: Arrays, like other Python sequences can also be sliced >>> a = mp-arange(10) array({0, 1, 2, 3, 4, 5, 6, 7, 8 91) >>> al2:9:3} # [start :end-step) array([2, 5, 8)) Note that the last index is not included! >>> alsa] array([0, 1, 2, 31) All three slice components are not required: by default, start is 0, end is the last and step is 1 >>> alta) array([t, 21) >>> als2) array([0, 2, 4, 8, 81) >>> al3:] array({3, 4, 5, 6, 7, 8, 91) A small illustrated summary of NumPy indexing and slicing, 4.1, The NumPy array object a7Scipy lecture notes, Edition 2022.1 >>> a[®, 3:5] array((3, 4]) pep al4:, 4:] array([[44, 55], (54, 55]]) >e> al, 2] a([2, 12, 22, 32, 42, 521) pep a[2::2, 1:2] array({[26, 22, 24], (40, 42, 44]]) You can also combine assign: ant and slicing: >>> a = mp-arange(10) >>> als:] = 10 array([0, 1, 2, 3, 4 10, 10, 10, 10, 10]) >>> b = mp. arange(S) pop alB:] = Est] array([0, 1, 2, 3, 4, 4, 3, 2, 1, oD Exercise: Indexing and slicing + Try the different flavours of slicing, using start, end and step: starting from a linspace, try to obtain odd numbers counting backwards, and even numbers counting forwards, + Reproduce the slices in the diagram above. You may use the following expression to create the >>> mp.arange(6) + ap.arange(0, 51, 10)[:, np.newaxia) array([l0, 1, 2, 3, 4, 51, (10, 11, 12, 13, 14, 15), (20, 21, 22, 23, 24, 25), (30, 31, 32, 33, 34, 36), (ao, 41, 42, 43, 44, 457, (50, 51, 52, 53, 54, 58]]) Exercise: Array creation Create the following arrays (with correct data types): ta, 1, 1a, tt, 1, a0, (1,1, 21, (1,6, 1 1 (fo. 0., 0,0, 0.1, (2, 0.,0.,.0., 0.1, (0., 3.,0.,0., 0.1, (0. 0, 4, 0., 0.1, (0.,0.,0.,5., 0.1, (0.,0.,0.,0., 6.1] 4.1, The NumPy array object 58Scipy lecture notes, Edition 2022.1 Par on course: 3 statements for each Hint: Individual array elements can be accessed similarly to a list, e.g. alt] or alt, 2) Hint; Examine the docstring for diag, Exercise: Tiling for array creation Skim through the documentation for np. tile, and use this funetion to construct the array: (4, 3,43, 4 3), (2. 1,2, 4,2, 41, (4.3, 4,3, 4, 31, (1, 2,4, 2, 111 4.1.6 Copies and views A slicing operation creates a view on the original array, which is just a way of accessing array data. Thus the original array is not copied in memory. You can use np-may_share_memory() to check if two arrays share the same memory block. Note however, that this uses heuristics and may give you false positives, When modifying the view, the original array is modified as well. 55> a= mp.arange(i0) array([0, 1, 2, 3, 4, po b= als2] >>> b array({0, 2, 4, 6, 81) may_share nenory(a, b) 67,8, 91) >>> pio] = 12 array([12, 2, 4 pra #0) array(12, 1, 2, 3, 4, 5, >>> a= mp. arange(10) doe = ali:2] copy # force « copy >o> elo] = 12 array([0, 1, 2, 3, 4, 5, 8, 7, 8, SD) >>> mp.may_share_menory(a, ¢) False This behavior can be surprising at first sight... but it allows to save both memory and time Worked example: Prime numb: 4.1, The NumPy array object 59Scipy lecture notes, Edition 2022.1 > 01234567 8 910111213 3 RECORD RCERKEL BEE --. 4 RECCEC ED REDE «. VN Compute prime numbers in 0-99, with a sieve + Construct a shape (100,) boolean array is_prime, filled with True i the beginning: >>> is_prime = np.ones((100,), dtype=bool) + Cross out 0 and 1 which are not primes: >>> 1s_primel:2] = 0 + For each integer j starting from 2, cross out its higher multiples: » >> Wmax = int (ap. sqrt (len is_prime) — >>> for j im range(2, Nmax ~ 1) is prime(2+j!:3] = False + Skim through help(np.nonzero), and print the prime numbers + Follow-up! = Move the above code into a script file named prire_sieve.py — Run it to check it works ‘the optinization suggested in the sieve of Eratosthenes: 1. Skip j which are already known to not be primes 2. The first number to cross out is j2 4.1.7 Fancy indexing Tip: NumPy arrays can be indexed with slices, but also with boolean or integer arrays (masks). This method is called fancy indexing. It creates copies not views. Using boolean masks 5>> mp randon. seed(3) >>> a= np.random.randint(0, 21, 19) array([10, 3, 8, 0, 19, 10, 11, 9, 10, 6, 0, 20, 12, 7, 141) >> 3-20) array([False, True, Felse, True, Felse, False, False, True, False (conties om HER DaE) 4.1, The NumPy array object 60Scipy lecture notes, Edition 2022.1 {continued from previous page) False, True, False, False) °) almask] # or, 2 Lo% >>> extract from_a # extract a sub-array with the mask array([ 3, 0, 8, 8, 0, 121) de> mask = (a3 >>> extract from_a Indexing with a mask can be very useful to assign a new value to a sub-array: de> ala 3 s= 0) = = array([10, “1, 8, 1, 19, 10, 11, -1, ot, -1, 20, -1, 7, 14) Indexing with an array of integers 55> a= mp.arange(0, 100, 20) array({ 0, 10, 20, 30, 40, 50, 60, 70, 80, 901) Indexing can be done with an array of integers, where the same index is repeated several time: pe» alla, 3,2, 4, 21) # note: (2, 3, 2, 4, 2] 1s @ Python list array([20, 30, 20, 40, 20]) New values can be assigned with this kind of indexing: >> alls, 71 array({ 0, 10, 20, 30, 40, 60, 60, -100, 80, -100}) 100 Tip: When a new array is created by indexing with an array of integers, the new array has the same shape as the array of integers: 55> a= mp arange(10) de> idx = mp.array([I3, J, (8, 711) >>> idx. shape 2, 2 >>> alie) array (C13, 4, (9, 71) ‘The image below illustrates various fancy indexing applications >>> al(®,1,2,3,4), (1,2,3,4,5)] array({1, 12, 23, 34, 45]) >>> al3:, [0,2,51] array({[30, 32, 35], (40, 42, 45), [50, 52, 55]]) pop mask = np.array([1,0,1,0,0,1], dtype=bool) >>> almask, 2] array([2, 22, 521) 4.1, The NumPy array object 61Scipy lecture notes, Edition 2022.1 Exercise: Fancy indexing + Again, reproduce the fancy indexing shown in the diagram above, + Use fancy indexing on the left and array ereation on the right to assign values into an array, for instance by setting parts of the array in the diagram above to zero. 4.2 Numerical operations on arrays Section contents ise operations + Basic reductions + Broadcastins + Array shape manipulation + Sorting data + Summary 4.2.1 Elementwise operations Basic operations With scalars: >>> a = np.array(li, 2, 3, 4) ways, >>> Drea, array([2, 4, 8, 16]) 4,81) All arithmetic operates elementwise: >>> b= mp.ones(a) + 1 array([=1., 01, dey 2.1) p> asd array(I2., 4, 6, 8.J) >>> j = mp-arange(s) de> Beet = 5 array(L 2, 3, 8, 13, 28)) ‘These operations ave of course much faster than if you did them in pure python: 35> a = np. arange(10000) do> Yeimest a+ 10000 loops, best of 3: 24.3 us per loop >>> 1 = range(10000) >>> eimai [iel for 4 in 2] 1000 loops, best of 3: 861 us per loop Warning: Array multiplication is not matrix multiplication: 4.2. Numerical opera 62Scipy lecture notes, Edition 2022.1 >>> ¢ = mp.onea((3, 3)) array, 1, 1d, Qo 4; 12 G4, 19D Note: Matrix multiplication: >>> €.dor(e) array([(3., 3., 3.1, sven stent B., 3, 31D xercise: Elementwise operations + Try simple arithmetic elementwise operations: add even elements with odd elements + Time them against their pure python counterparts using /timest, + Generate: = (2000, 2001, 200! 243, 2804) nape 2 ep - 5 Other operations Comparisons: de> a= mp.array([1, 2, 3, 41) de> b= mp.array([4, 2, 2, 41) poo ase array([Faise, True, False, True) p> ard array([False, False, True, False]) Tip: Array-wise comparisons: >>> a= mp.array([1, 2, 3, 4)) de> b= mp.array((a, 2, 2, 4) de> e = mp.array((1, 2, 3, 41) >>> ay_equal(a, b) False >>> mp.array_equal(a, ¢) True Logical operations: >>> a= mp.array(L1, 1, >>> b= mp.array([I, 0, >>> mp-logical_er(a, ) array({ True, True, True, False]) >>> mp-logical_ana(e, b) array({ True, False, False, False) ‘Transcendental functions: 4.2. Numerical opera 63Scipy lecture notes, Edition 2022.1 >>> a = mp.arange(S) >>> mp.ein(a) array(L 0 + 0.84247098, 0.90929743, 0.14112001, -0.7568025 1) >>> mp-Leg(a) array int, 0 |» 0.89814718, 1,09861229, 1.38629436)) >>> mp. exp(a) array(l 1 + 2.71828183, 7.890561 , 20,08853682, 54.59815003)) Shape mismatches 55> a = mp.arange() >>> a+ mp-array([1, 21) Traceback (most recent call last) File "", Line 2, in ValueError: operands could not be broadcast together with shapes (4) (2) Broadcasting? Well return to that later ‘Transposition >>> a= mp.criu(mp.ones((3, 3)), 1) # see help(np. tris) array(C[0., 1, 1.1, fo, 0. 10, fo, 0, ol >> at array({[0., 0., 0-1, Ley : fy a, 01D Note: The transposition is a view The transpose returns a view of the original array >>> a = up.arange(9).reshape(3, 5) >>> a.Tlo, 2] = 999 >> at array([{ 9, 3, 999], Cu oa om [2 5, sid array([l 0, 4, 2, C3 4 8h (999, 7, 81]) Note: Linear algebra The sub-module nunpy. linalg implements basic linear algebra, such as solving linear syster value decomposition, etc. However, it is not guaranteed to be compiled using eflicien routines, and singul we recommend the use of scipy.1inalg, as detailed in section Linear algebra operations: seipy.tinaly Exercise other operations + Look at the help for np-allclose, When might this be useful? + Look at the help for np.triu and ap.tril 4.2. Numerical opera 64Scipy lecture notes, Edition 2022.1 4.2.2. Basic reductions Computing sums >o> x= mp.array((1, 2, 3, 40) >>> ap. sunx) 10 >>> x. sun() 10 Sum by rows and by columns: axis 0 >> x= aparray([l, 0, 2, 2D array( n, (2, 21) >>> x, sun (axie-0) array({3, 3)) >o> xfs, O)-sunQ, xfs, 1) sun @, 3) >>> x.sun(axi array((2, 41) >>> xf0, :J-gumQ, x[t, :J-sumQ 2,4 # columns (frst dimension) 1) # rows (second dimension) Tip: Same idea in higher dimensions: >>> x = np.random.rand(2, 2, 2) >>> x.sun(axis=2)[0, 1] 1.764. >>> xf0, 1, :].sumQ) 1.14764 Other reductions — works the same way (and take axis=) Extrema: >>> x = mp.array([1, 3, 20) >>> x.min() 1 >>> x.mex() (conties om He Dae) 4.2. Numerical operations on arrays cSScipy lecture notes, Edition 2022.1 3 de> xargninQ) # index of mintmum ° >>> x.argnax() # indes of maztmun 1 Logical operations: ali((irue, True, False]) Trve, False]) {continued from previous page) Note: Can be used for array comparisons >>> a= mp.zeres((100, 100)) >>> ap.any(a '= 0) False >>> mp.all(a == a) True de> a= mp-array((1, 2, 3, 23) >>> b de> c= mp.array((5, 4, 4, 5) d>> (la @ b) tb & )) allO true Statistic de> x > mp array((1, 2, 3, 10) >>> y= mp.array(([Ul, 2, 31, (5, 5, 11) >>> x.meanQ) 1.78 >>> mp.median(x) 15 >>> mp.median(y, axis: array([2., 5.1) b> x. sta) ¥ full population standard dev 0.829156:9758884995 and many more (best to leam as you go). Exercise: Reductions + Givon there is a sum, what other function might + What is the difference between sum and cursum? you expect to see? Worked Example: diffusio using a random walk algorithm 4.2. Numerical opera 66Scipy lecture notes, Edition 2022.1 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 Fight jumps? We are going to sinilate 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: shuffled jumps Position : cumulated jumps sum time time 1j1}-1]1}-2)1 1) 2|1|2{1/2 -1)1]1}2}1/1 -1/0/}1/0/1/0 » o a a S S a a >> m_ateries = 1000 # nunber of vetkers >>> tina = 200 # tine during which we follow the walker ‘We randomly choose all the steps 1 or -1 of the walk: >>> t = mp arange(t_nax) >>> steps = 2 * mp.random.randint (0, 1+ 1, (n_storie value is exclusive >>> mp.unique(ateps) # Verification: all steps ere 1 or ~ array((-1, 1) ax) - 1 # +1 because the highs We build the walks by summing steps along the time: >> positions = np.cunsun(steps, axis=1) # asic = >>> aq distance = positions*#2 ‘We get the mean in the axis of the st inate aT Plot the results: >>> pt figure(figsize-(4, 3))

de> plt.plot(t, np.sqrt(mean_sq_distance), ‘g.', t, np.sqrt(t), 'y-') (enatplotlib.LinesLine2D object at ...>, ) de> pt. xlabel (r'S¢ Toxt(. $8") de> pt. ylabel(r"S\sqre(\langle (\delta x)°2 \rangle} ‘Text (..."$\\sqrt(\\langle (\\delta x)-2 \\rangle)$") >>> plt tight layout() # provide sufficient space for labels ) > 4.2. Numerical operations on arrays o7Scipy lecture notes, Edition 2022.1 100150 t We find a well-known result in physics: the RMS distance grows as the square root of the time! 4.2.3 Broadcasting + Basic operations on numpy arrays (addition, etc.) are elementwise + This works on arrays of the sane size Nevertheless, It’s also possible to do operations on arrays of different sizes if NumPy can transfon ‘these arrays so that they all have the same size: this conversion is called broadcasting, The image below gives an example of broadcasting: ofofo [42 ofofo (4 [2] soltofto], , ofaf2|) _ |rofrofio)| , [olajz 20 [20 |20 ofaf2 20|20|20 oft|2 \ 30/30 [30 ef 30| 3030 oft|2 N o[ojo o[+[2 ofolo ORE alalz ‘tol s9 [901 . + = [ropa faz 20|20| 20 20/21 [22 30]30[30 3031 | 32 ° o(i{2 ° oft(z]) vi 10 7 = | 7 7 20 20 30 30 Let’s verify: 4.2. Numerical operations on arrays 68Scipy lecture notes, Edition 2022.1 >>> a= mp.cile(mp.arange(0, 40, 10), (3, 1)).7 array({L 0, 0, ol, 10, 10, 10), (20, 20, 20), [30, 30, 3011) de> b= mp array([0, 1, 2]) pat array([l.0, 1, 21 Lio, 11, 121 (20, 21, 221, [30, 31, 3211) We have already used broadcasting without knowing it! >>> a= mp ones((4, 5)) >>> alo] = 2 # ve assign an array of dimension 0 to an array of dimension 1 array({[2., 2., 2 241 oe 41, re 41, a 1a A useful trick: >>> a= mp.arange(0, 40, 10) >>> a. shape a >>> a= al:, mpnevaxis] # adds a new azis ~> 20 array >>> a.shape @p array({{ 0), array([l 0, 1, 21, 10, 11, 121, (20, 21, 221, (30, 31, 3211) ‘Tip: Broadcasting seems a bit magical, but it is actually quite natural to use it when we want to solve a problem whose output data is an array with move dimensions than input data, ‘Worked Example: Broadcasting Let’s construct an array of distances (in miles) between cities of Route 66: Chicago, Springfield, Saint-Louis, Tulsa, Oklahoma City, Amarillo, Santa Fe, Albuquerque, Flagstaff and Los Angeles. >> milleposts = ap.array((0, 198, 303, 736, 871, 1175, 1475, 1844, 1913, 2448)) de> distance array ~ mp.abs(aileposts ~ mileposts[:, np nevaxis]) >>> distance array array([l 0, 198, 303, 736, 871, 1275, 1475, 1844, 2s], L198, 0, 105, 538, 673, 977, 1277, 1345, 2250], [ 303, 105, 0, 433, 568, 872, 1172, 1241, 2610, 2145], [ 738, 538, 433, 0, 135, 439, 739, 08, 17121, Lan, 568, 135, 0, 604, 873, 1871, (ur, 37254395308 TET TT. Tse, oo 5 [is44, 1348, 1241, 808, 673, 369, 69, 0,369, _ soa] 4.2. Nubieittaltoperatons dA7derayl2, 738, 438, 369, 0, S36), Gi (244s, 2250, 2145, 1712, 1577, 1273, 973, 904, 535, O11) 300,369, 1273)Scipy lecture notes, Edition 2022.1 A lot of grid-based or network-based problems can also use broadcasting. For instance, if we want to compute the distance from the origin of points on a 5x5 grid, we can do >>> x, y = mp-arange(S), np-arange(s)[:, mp.nevaxis] >o> distance ~ np.sqrt(x “+ 2 + y + 2) >>> distance array (L(0. al o 8 o & 4 1, ree | 441421366, 2.23606798, 3.16227766, 4. 12310563! i } 2.23606788, 2.82842712, 3.60555128, 4.47213595; b& 3.60555128, 4.464059, 5. le 447213595, 5. + 5.85685425; Or in color: >>> ple peoler (distance) >o> ple. colerbar() Remark : the nunpy-ogrid() function allows to directly create vectors x and y of the previous example, with two “significant dimensions”: Do x, y = mpiogria[0:5, 0:5) >> xy (array (C10) , a, a, BI, (41), array(C0o, 1, 2, 3, 41D) (conties om He Dae) 4.2. Numerical opera 70Scipy lecture notes, Edition 2022.1 {continued from previous page) >>> x.shape, y shape: S, 1), G89) >> distance = np.sqrt(x “+ 2 + y *¥ 2) So, np. ogrid is very useful as soon as we have to handle computations on a grid. On the other hand, np-mgrid directly provides matrices full of indices for cases where we can’t (or don't want to) benefit from broadcasting, >>> x, y= mp.mgrid[0r4, 0:4) array({{0, 0, 0, Oo), ena, Iz, 2, 2, 21, (3, 3,3, 311) > array({00, 1, 2, 31, [o, 1, 2, 31, (0, 1, 2, 31, fo, 1, 2, 311) See also: Broadcasting: discussion of broadcasting, in the Advanced NurPy chapter. 4.2.4 Array shape manipula Flattening de> a= mp.array((ll, 2, 3], (@, 5, 6l]) >>> acravel) array([t, 2, 3, 4, 5, 61) >> aT array({Lt, 4, isi, (3, 811) de> a.Turavel() array((t, 4, 2, 8. 3, 6)) Higher dimensions: last dimensions ravel out “first” Reshaping ‘The inverse operation to flattening: 55> a. shape @, 3) >>> b= a-ravel >>> b= boreshape((2, 3)) >>> b array({{1, 2, 31, [4 5, 61) Or, >>> areshape((2, -1)) # unspecified (-1) value as inferred array(EO3, 2, 3], la, 5, 611) 7 4.2. Numerical operaScipy lecture notes, Edition 2022.1 Warning: ndarray.reshape may return a view (cf help(np. reshape))), or copy Tip: >>> BIO, 0) = 98 array(((98, 2, 31, la, 5, 611d Beware: reshape may also return a copy! >>> a= mp.zeres((S, 2)) >>> b= aT reshape(s+2) >>> bio] = 8 array({(0., 0] f0., 0.3 o., 0.1 2 ‘To understand this you need to learn more about the memory layout of a numpy array. ‘Adding a dimension Indexing with the ap.nevaxis object allows us to add an axis to an array (you have seen this alzeady above in the broadcasting section): >>> z= mp array([1, 2, 30) array((t, 2, 3)) de> ali, ap nevaxis] array (CLI, (2 (31) >>> almp.nevaxis, :1 array(C[3, 2, 31)) Dimension shuffling de> a= mp. arange(493+2).reshape(4, 3, 2) >>> a. shape (, 3, 2) >>> alo, 2, 1] 5 >>> b= actranspose(i, 2, 0) >>> shape «3, 2, 4) >>> bI2, 1, 0] 5 Also ereates a view: >>> BL, , ol >>> alo, 2, 1] a 2 4.2. Numerical operaScipy lecture notes, Edition 2022.1 Resizing Size of an array can be changed with ndarray resize: np. arange(4) ssize((8,)) array({0, 1, 2, 3, 0, 0, 0, 01) However, it must not be referred to somewhere else: b> baa >>> a.resize((é,)) ‘Traceback (ost recent atl last) File "", line 2, in ValueError: cannot resize an array that has been referenced or is referencing another array in this way. Use the resize function Exercise: Shape manipulations + Look at the docstring for reshape, especially the notes section which has some more information about copies and views. Use flatten as an alternative to rave). What is the difference? (Hint: check which one returns a view and which a copy) + Experiment with transpose for dimension shuflling 4.2.5 Sorting data Sorting along an axis: de> a= mp.array(([4, 3, 5], 1, 2, 10) >>> b= mp.sort (a, axis >>> array({[3, 4, 51, fa, ny Note: Sorts each row separately! In-place sort: de> a. gort (axisn1) array({(3, 4, 5] ty, a, 21 Sorting with fancy indexing: >o> a= mp.array((, 3, 1, 22) >>> j= mp-argzort (a) >> 5 array([2, 3, 1, 01) >>> all array([t, 2, 3, 41) Finding minima and maxima, 4.2. Numerical opera 73Scipy lecture notes, Edition 2022.1 >>> a= mp.array([4, 3, >>> j.max = np-argnax(a) >>> jumin = np.argnin(a) 2D Exercise: Sor + Thy both implace and out-of-place sorting. + Thy creating arrays with different dtypes and sorting them. se all or array equal to check the results Look at np. random. shuffle for a way to create sortable input quicker Combine ravel, sort and reshape. Look at the axis keyword for sort and rewrite the previous exercise 4.2.6 Summary What do yc need to know to get started? + Know how to create arrays : array, arange, ones, zeros. + Know the shape of the array with array shape, then use slicing to obtain different views of the anay: array(::21, ete. Adjust the shape of the array using reshape or flatten it with ravel. + Obtain a subset of the elements of an array and/or modify their values with masks de> ala <0] =0 + Know miscellaneous operations on arrays, such as finding the mean or max (array.max(), array snean()). No need to retain everything, but have the reflex to search in the documentation (online does, helpQ, leokfer())! + For advanced use: master the indexing with arrays of integers, as well as broadcasting. Know more NumPy functions to handle various array operat Quick read If you want to do a first quick pass through the Scipy lectures to learn the ecosystem, you can directly skip to the next chapter: Matplotlib: pl The remainder of this chapter is not necessary to follow the rest of the intro part. But be sure to come back and finish this chapter, as well as to do some more exercices 4.3 More elaborate arrays Section contents + More ata types + Structured data types + maskedarray: dealing with (propagation of) missing data 4.3. More elaborate arrays 74