<jupyter_start><jupyter_text>NumPy -  multidimensional data arrays ... with some Matplotlib`<NAME>`Based on the work of- <NAME> (MNE) and- <NAME> (<EMAIL>) http://dml.riken.jp/~rob/<jupyter_code># what is this line all about?!? Answer later
%matplotlib inline<jupyter_output><empty_output><jupyter_text>Introduction The `numpy` package (module) is used in almost all numerical computation using Python. It is a package that provide high-performance vector, matrix and higher-dimensional data structures for Python. It is implemented in C and Fortran so when calculations are vectorized (formulated with vectors and matrices), performance is very good. To use `numpy` first import the module:<jupyter_code>import numpy as np # np is the most common import name for numpy<jupyter_output><empty_output><jupyter_text>In the `numpy` package the terminology used for vectors, matrices and higher-dimensional data sets is *array*.   Creating `numpy` arrays There are a number of ways to initialize new numpy arrays, for example from* a Python list or tuples* using functions that are dedicated to generating numpy arrays, such as `arange`, `linspace`, etc.* reading data from files  From lists For example, to create new vector and matrix arrays from Python lists we can use the `numpy.array` function.<jupyter_code># a vector: the argument to the array function is a Python list
v = np.array([1, 0, 3, 4])
# print v, and the type of v
print(v)
type(v)
# a matrix (or better a 2d array): the argument to the array function is a nested Python list
M = np.array([[1, 0], [3, 4]])
# print M
print(M)<jupyter_output>[[1 0]
 [3 4]]<jupyter_text>The `v` and `M` objects are both of the type `ndarray` that the `numpy` module provides.<jupyter_code># show the type of v and of M
type(M)<jupyter_output><empty_output><jupyter_text>The difference between the `v` and `M` arrays is only their shapes. We can get information about the shape of an array by using the `ndarray.shape` property.<jupyter_code># get the shape of v
v.shape
# note the type of v.shape
type(v.shape)
# get the shape of M
print(M.shape)<jupyter_output>(2, 2)<jupyter_text>The number of elements in the array is available through the `ndarray.size` property:<jupyter_code># get the size -- be careful matlab size is shape with NumPy!
print(v.size)
print(M.size)<jupyter_output>4
4<jupyter_text>So far the `numpy.ndarray` looks awefully much like a Python list (or nested list).**Why not simply use Python lists for computations instead of creating a new array type?** There are several reasons:* Python lists are very general:  - They can contain any kind of object  - They are dynamically typed  - They do not support mathematical functions such as matrix and dot multiplications, etc. NumPy array are:* **statically typed** and **homogeneous**. The type of the elements is determined when array is created.* Memory efficient.* Allow fast implementation of mathematical functions in compiled code.Using the `dtype` (data type) property of an `ndarray`, we can see what type the data of an array has:<jupyter_code>print(M.dtype) # show M's dtype<jupyter_output>int64<jupyter_text>We get an error if we try to assign a value of the wrong type to an element in a numpy array:<jupyter_code>M[0, 0] = 'bad'  # try to assign "hello" to the upper left element
M[0, 0] = 1.5 # try to assign a float to the upper left element
print(M) # print the result<jupyter_output>[[1 0]
 [3 4]]<jupyter_text>If we want, we can explicitly define the type of the array data when we create it, using the `dtype` keyword argument:<jupyter_code>M = np.array([[1, 0], [3, 4]], np.complex) # create M, a complex array with shape (2, 2)
M[0, 0] = 1.5
print(M) # print M<jupyter_output>[[1.5+0.j 0. +0.j]
 [3. +0.j 4. +0.j]]<jupyter_text>Common types that can be used with `dtype`:* `int`* `float`* `complex`* `bool`We can also explicitly define the bit size of the data types, for example: `int64`, `int16`, `float128`, `complex128`.  Using array-generating functions For larger arrays it is inpractical to initialize the data manually, using explicit python lists. Instead we can use one of the many functions in `numpy` that generates arrays of different forms. Some of the more common are:  arange<jupyter_code>x = np.arange(0, 10)  # create x from 0, 1, ..., 9
print(x)
x = np.arange(-1, 1, 0.1)  # create x from -1 to 1 in steps of 0.1
print(x)<jupyter_output>[-1.00000000e+00 -9.00000000e-01 -8.00000000e-01 -7.00000000e-01
 -6.00000000e-01 -5.00000000e-01 -4.00000000e-01 -3.00000000e-01
 -2.00000000e-01 -1.00000000e-01 -2.22044605e-16  1.00000000e-01
  2.00000000e-01  3.00000000e-01  4.00000000e-01  5.00000000e-01
  6.00000000e-01  7.00000000e-01  8.00000000e-01  9.00000000e-01]<jupyter_text>linspace and logspace<jupyter_code>x = np.linspace(0, 1, 10)  # using linspace, both end points ARE included
print(x)
x = np.linspace(0, 10, 25)  # create array from 0 to 10, 25 elements
print(x)
# create log-spaced data, 0 to 10, 10 elements
print(np.logspace(0, 10, 10))<jupyter_output>[1.00000000e+00 1.29154967e+01 1.66810054e+02 2.15443469e+03
 2.78255940e+04 3.59381366e+05 4.64158883e+06 5.99484250e+07
 7.74263683e+08 1.00000000e+10]<jupyter_text>First demo of MatplotlibFirst, import:<jupyter_code>import matplotlib.pyplot as plt # import matplotlib is most often done using "plt" as an alias
x = np.arange(0, 10)  # x data from 0 to 9
y = np.linspace(0, 1, 10) # y data as 0 to 1, 10 elements
plt.plot(x, y, 'r-')  # plot the data as a red line
plt.show()  # show the plot<jupyter_output><empty_output><jupyter_text>random data<jupyter_code>from numpy import random # import random submodule of NumPy
x = random.rand(5, 5) # shape (5, 5) uniform random numbers in [0,1]
print(x)
x = random.randn(5, 5)  # standard normal distributed random numbers
print(x)<jupyter_output>[[ 0.18057247 -0.92300503 -0.54869637  1.138543   -0.06970905]
 [ 0.49238966  0.7861138   1.48629807  2.21248713 -0.18132759]
 [ 0.57581772 -1.1625285  -0.05458411 -0.43917363  0.28151315]
 [ 0.34397356  0.25352392 -0.78481538 -0.8754882   1.13391234]
 [-0.71792639 -0.40699086 -0.06653698 -0.15145506  0.84504088]]<jupyter_text>More matplotlib<jupyter_code>x = random.randn(100000)  # 100000 normally distributed numbers
plt.hist(x)  # matplotlib hist(ogram)
plt.show()  # show it
x = random.randn(10, 10) # shape (10, 10) random data
plt.imshow(x, interpolation='nearest', cmap='viridis')  # imshow with nearest interpolation, cmap viridis
plt.colorbar() # add a colorbar
plt.show()  # show it<jupyter_output><empty_output><jupyter_text>diag<jupyter_code>y = np.diag(np.diag([1, 3, 5]))  # create a diagonal matrix
print(y)
y = np.diag([1, 3, 5], k=-1)  # diagonal with offset (k) from the main diagonal 
print(y)<jupyter_output>[[0 0 0 0]
 [1 0 0 0]
 [0 3 0 0]
 [0 0 5 0]]<jupyter_text>zeros and ones<jupyter_code>x = np.zeros((3, 3)) # create a an array of shape (3, 3) of zeros
x[:] = 5
print(x)
x = np.ones((3, 3))  # similar for ones
print(x)
x = np.full((3, 3), 5)  # and for any single value (full)
print(x)<jupyter_output>[[5 5 5]
 [5 5 5]
 [5 5 5]]<jupyter_text>File I/O Using `numpy.savetxt` we can store a Numpy array to a file in CSV format:<jupyter_code>M = random.randn(3, 3)  # M as (3, 3) rand
np.savetxt('temp.txt', M)  # savetxt writes these to a text file with reasonable defaults
!cat temp.txt # use a system command (!cat) to check
np.savetxt('temp.txt', M, fmt='%0.1f', delimiter='\t') # change fmt and delimiter
!cat temp.txt # look at the result with !cat
# load with loadtxt, minding the delimiter
y = np.loadtxt('temp.txt')
print(y)<jupyter_output>[[ 2.3  0.6 -0.3]
 [-0.9 -0.9  2.1]
 [-0.8  0.3  0.3]]<jupyter_text>Numpy's native file format Useful when storing and reading back numpy array data. Use the functions `numpy.save` and `numpy.load`:<jupyter_code>np.save('temp.npy', M) # np.save writes in binary format
!file temp.npy  # inspect using !file
M2 = np.load('temp.npy')  # np.load loads it again
print(M, M2)<jupyter_output>[[ 2.3344224   0.5966002  -0.32530528]
 [-0.87824807 -0.9080117   2.09421343]
 [-0.76706903  0.26801627  0.28815671]] [[ 2.3344224   0.5966002  -0.32530528]
 [-0.87824807 -0.9080117   2.09421343]
 [-0.76706903  0.26801627  0.28815671]]<jupyter_text>More properties of the numpy arrays<jupyter_code>print(M.dtype) # dtype
print(M.itemsize) # itemsize: bytes per element
print(M.nbytes)  # nbytes: number of bytes
print(M.size * M.itemsize)  # ... which is the size * itemsize
print(M.ndim, len(M.shape))  # ndim: number of dimensions<jupyter_output>2 2<jupyter_text>Manipulating arrays  Indexing We can index elements in an array using the square bracket and indices:<jupyter_code>v = np.array([1, 0, 3, 4])
print(v.shape)  # v is a vector, and has only one dimension, taking the first index
print(v[0])
print(M[0, 0])  # M is a matrix, or a 2 dimensional array, taking two indices<jupyter_output>2.3344224007974503<jupyter_text>If we omit an index of a multidimensional array it returns the whole row (or, in general, a N-1 dimensional array)<jupyter_code>print(M) # print all of M
print(M[0])  # print row 1 of M<jupyter_output>[ 2.3344224   0.5966002  -0.32530528]<jupyter_text>The same thing can be achieved with using `:` instead of an index:<jupyter_code>print(M[0, :])  # row 1
print(M[:, 0]) # column 1<jupyter_output>[ 2.3344224  -0.87824807 -0.76706903]<jupyter_text>We can assign new values to elements in an array using indexing:<jupyter_code>M[0, 0] = 5. # assign to upper left of M
print(M) # print M
# also works for rows and columns
M[1] = 0.5 # assign to row 1 of M
M[:, 0] = 2 # assign to column 2 of M
print(M)
# print M<jupyter_output><empty_output><jupyter_text>Index slicing Index slicing is the technical name for the syntax `M[lower:upper:step]` to extract part of an array:<jupyter_code>A = np.array([1, 2, 3, 4, 5])
print(A)
print(A.shape)
A[1:3]  # get an array slice from index 1 to 3 (excluding endpoint!)<jupyter_output><empty_output><jupyter_text>Array slices are *mutable* because they create a *view* of the array, rather than a *copy*: if they are assigned a new value the original array from which the slice was extracted is modified:<jupyter_code>B = A[1:3]  # assign B part of A
B[0] = 100  # modify all of B
print(A) # show that it affects A<jupyter_output>[  1 100   3   4   5]<jupyter_text>We can omit any of the three parameters in `M[lower:upper:step]`:<jupyter_code>A[::2] # with : or :: t hen lower, upper, step all take the default values
# step is 2, lower and upper defaults to the beginning and end of the array
A[:3]  # first three elements
A[3:] # elements from index 3 onward
B = np.zeros((4, 4))
B[1] = 1
print(B)
print(B[1::2, :])<jupyter_output>[[0. 0. 0. 0.]
 [1. 1. 1. 1.]
 [0. 0. 0. 0.]
 [0. 0. 0. 0.]]
[[1. 1. 1. 1.]
 [0. 0. 0. 0.]]<jupyter_text>Negative indices count from the end of the array (positive index from the begining):<jupyter_code>A = np.arange(1, 6) # create A as a vector from 1 to 5
print(A)
A[-1]  # get the last (-1) element in the array
A[-3:]  # get the last three elements<jupyter_output><empty_output><jupyter_text>Index slicing works exactly the same way for multidimensional arrays:<jupyter_code>x = np.array([[10 * row + col for col in range(1, 5)] for row in range(1, 5)])  # create 4x4 with tens digit as row, ones digit as col index
print(x)
# show middle 2x2
print(x[1:3, 1:3])
# stride by 2 in each dimension
print(x[::2, ::2])<jupyter_output>[[11 13]
 [31 33]]<jupyter_text>Fancy indexing Fancy indexing is the name for when an array or list is used in-place of an index:<jupyter_code>A = random.randn(4, 4)
print(A[[0, 1, 2]])
# Get 3 elements, these rows matched these with
# cols 0, 1, last (recall -1 means the last element)
A[np.ix_([0, 3], [0, 3])]
# to actually do row x col combinations, you can np.ix_ them<jupyter_output><empty_output><jupyter_text>We can also index masks: If the index mask is an Numpy array of with data type `bool`, then an element is selected (True) or not (False) depending on the value of the index mask at the position each element:<jupyter_code>B = np.arange(0, 5)  # create B as the integers from 0 to 4
mask = B % 2 == 0  # get the indices of even elements using mod 2 
print(mask)  # manually create a 5-element mask using bools, use it
B[[True, False, True, True, True]]
B[np.array([1, 0, 1, 1, 1], bool)]  # same thing with np.array(..., dtype=bool) casting of ints<jupyter_output><empty_output><jupyter_text>This feature is very useful to conditionally select elements from an array, using for example comparison operators:<jupyter_code>x = np.linspace(0, 1, 10)  # create an array x using linspace
mask = x > 0.5  # create a mask for elements between 5 and 7
x[mask] # apply this mask to x<jupyter_output><empty_output><jupyter_text>Functions for extracting data from arrays and creating arrays  where The index mask can be converted to position index using the `where` function<jupyter_code>np.where(mask)[0] # get indices from mask
x[np.where(mask)[0]]  # index x using indices: equivalent to the fancy indexing x[mask]<jupyter_output><empty_output><jupyter_text>diag With the diag function we can also extract the diagonal and subdiagonals of an array:<jupyter_code>print(np.diag(A))  # get diagonal elements of A
print(np.diag(A, -1))  # get diagonal elements offset left by one<jupyter_output>[ 0.68676875  0.17069703 -1.41924581]<jupyter_text>Array oriented algebra 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.  Scalar-array operations We can use the usual arithmetic operators to multiply, add, subtract, and divide arrays with scalar numbers.<jupyter_code>v1 = np.arange(1, A.shape[0] + 1)  # create v1 as integers 0 to A.shape[0] - 1
print(v1)
v1 * 2# multiply by two
v1 + 2# add two<jupyter_output><empty_output><jupyter_text>Element-wise array-array operations When we add, subtract, multiply and divide arrays with each other, the default behaviour is **element-wise** operations:<jupyter_code>A * A # multiply A by itself: element-wise multiplication
# this is NOT a matrix / dot product!)
v1 * v1 # multiply v1 by itself<jupyter_output><empty_output><jupyter_text>If we multiply arrays with compatible shapes, we get an element-wise multiplication of each row:<jupyter_code>print(A.shape, v1[:, np.newaxis].shape)  # look at the shape of A and v1
print(A * v1[:, np.newaxis]) # tile v1 to match the shape of A 
# use np.newaxis to make v1 have two dimensions
# print shapes of A and v1
# multiply A by v1 using tile
# multiple A by v1 using np.newaxis
# do an outer product of v1 with itself using np.outer
# show that this is the same as using np.newaxis
# show that v1.reshape can be used to add axes<jupyter_output><empty_output><jupyter_text>Matrix algebra What about matrix mutiplication? There are two ways. We can either use the `dot` function, which applies a matrix-matrix, matrix-vector, or inner vector multiplication to its two arguments:<jupyter_code>np.dot(A, A) # dot A with itself
np.dot(A, v1)  # dot A with v1
# dot v1 with v1<jupyter_output><empty_output><jupyter_text>Data processing Often it is useful to store datasets in Numpy arrays. Numpy provides a number of functions to calculate statistics of datasets in arrays.<jupyter_code>A.shape # remind us of the shape of A<jupyter_output><empty_output><jupyter_text>mean<jupyter_code>A.mean()  # get the mean of A
A[:, -1].mean()  # get the mean of the last column of A
A.mean(axis=1)  # get the mean of all columns of A<jupyter_output><empty_output><jupyter_text>standard deviation, variance, min, max, peak-to-peak...<jupyter_code>np.ptp(A, axis=0)  # np.std, np.var, np.min, np.max, np.ptp
np.percentile(A, [10, 90])  # percentile<jupyter_output><empty_output><jupyter_text>sum, prod, and trace<jupyter_code>np.sum(np.arange(1, 11)) # create an array d of values 1 to 10
np.sum(A)  # sum up all elements
np.sum([1, 2, 5])

np.cumsum(np.arange(5))  # cummulative sum<jupyter_output><empty_output><jupyter_text>EXERCISE :Compute an approximation of $\pi$ using Wallis' formula (no loop just using Numpy): --><jupyter_code>my_pi = ...<jupyter_output><empty_output><jupyter_text>Calculations with higher-dimensional data When functions such as `min`, `max`, etc., is applied to a multidimensional arrays, it is sometimes useful to apply the calculation to the entire array, and sometimes only on a row or column basis. Using the `axis` argument we can specify how these functions should behave:<jupyter_code>A = random.randn(3, 3) # create a random array of shape (3, 3)
A.max()  # global max
A.max(axis=0)  # max in each column
A.max(axis=1)  # max in each row<jupyter_output><empty_output><jupyter_text>Many other functions and methods in the `array` and `matrix` classes accept the same (optional) `axis` keyword argument.  Reshaping, resizing and stacking arrays The shape of an Numpy array can be modified without copying the underlaying data, which makes it a fast operation even for large arrays.<jupyter_code>print(A)  # print A
np.reshape(A, (1, n * m))
np.reshape(A, (1, n * m)) # create array B that has shape (1, n * m)
B = np.reshape(A, (1, n * m)).copy() # modify B
B[0] = 10000 # print A, showing that it is also changed:
print(A)  # B is view of the same data!<jupyter_output>[[100. 100. 100.]
 [100. 100. 100.]
 [100. 100. 100.]]<jupyter_text>Adding a new dimension: newaxis With `newaxis`, we can insert new dimensions in an array, for example converting a vector to a column or row matrix:<jupyter_code># create array v of values 1, 2, 3
# show v's shape
# make a 1-column matrix of the vector v
# make a 1-row matrix of the vector v<jupyter_output><empty_output><jupyter_text>Stacking and repeating arrays Using function `repeat`, `tile`, `vstack`, `hstack`, and `concatenate` we can create larger vectors and matrices from smaller ones.**But in practice they are almost always possible to avoid**  tile and repeat<jupyter_code>a = np.arange(4)  # create an array a of shape (2, 2) of the values 1 through 4
np.repeat(a, 2)  # repeat each element twice
np.tile(a, 2)  # tile the matrix twice<jupyter_output><empty_output><jupyter_text>concatenate<jupyter_code>b = np.array([5, 6]) # make an array b of shape (2,): values 5 and 6
print(A[:2, :2].shape, b[:, np.newaxis].shape)
np.concatenate([A[:2, :2], b[:, np.newaxis]], axis=1)  # try to concatenate a and b
# fix b so that it's shape (1, 2)
# concatenate a and b arrays (default is axis=0)
# concatenate a and b.T along axis=1<jupyter_output><empty_output><jupyter_text>hstack and vstack<jupyter_code># vstack a and b
# hstack a and b.T<jupyter_output><empty_output><jupyter_text>Copy and "deep copy" To achieve high performance, assignments in Python usually do not copy the underlaying objects. This is important for example when objects are passed between functions, to avoid an excessive amount of memory copying when it is not necessary (techincal term: pass by reference).<jupyter_code># Create variable B that refers to the same array data as A 
A = random.randn(3, 3)
B = A.copy()
B[0, 0] = 100
print(A)
print(B)
B = A[0] # using single-index into A creates a view
B[0] = 1000
print(A)
# so modifying B[0] modifies A<jupyter_output><empty_output><jupyter_text>If we want to avoid this behavior, so that when we get a new completely independent object `B` copied from `A`, then we need to do a so-called "deep copy" using the function `copy`:<jupyter_code>B = A[0].copy()  # make B as a copy of A
B[0] = -100
print(A)
# now, if we modify B, A is not affected<jupyter_output><empty_output><jupyter_text>Iterating over array elements Generally, we want to avoid iterating over the elements of arrays whenever we can (at all costs). In interpreted languages 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:<jupyter_code>v = np.arange(1, 5)  # create v as values 1 through 4, iterate over it, print
for value in v:
    print(value)
M = random.randn(2, 2)  # reshape to make M of shape (2, 2), iterate over it twice
for row_idx, row in enumerate(M):
    print('row %d' % row_idx)
    for row_col in row:
        print(row_col)<jupyter_output>row 0
1.4381736362936857
-0.8831604668922024
row 1
-1.3146796360128634
-0.26781547244229326<jupyter_text>When we need to iterate over each element of an array and modify its elements, it is convenient to use the `enumerate` function to obtain both the element and its index in the `for` loop:<jupyter_code># iterate over the enumerate(M), squaring each element
# show that each element in M is now squared<jupyter_output><empty_output><jupyter_text>Using arrays in conditions When using arrays in conditions in for example `if` statements and other boolean expressions, one need to use one of `any` or `all`, which requires that any or all elements in the array evalutes to `True`:<jupyter_code>(M > 0).any(axis=1).all()
# use (M > ...).any() and print to show some logical branching
# use (M > ...).all(...) and print to show more branching<jupyter_output><empty_output><jupyter_text>Type casting Since Numpy arrays are *statically typed*, the type of an array does not change once created. But we can explicitly cast an array of some type to another using the `astype` functions (see also the similar `asarray` function). This always create a new array of new type:<jupyter_code>M.dtype  # show M's dtype
print(M)
M2 = np.round(M).astype(int)  # create M2 as int version
print(M2)
# show M2's dtype
M2.dtype
# create M3 as bool version<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.10914861321926314, 'fraction_numerical': 0.0635205372337979, 'mean_word_length': 4.745771670190275, 'pattern_counts': {'":': 0, '<': 198, '<?xml version=': 0, '>': 203, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 21, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '54', 'dir': 'jupyter-structured-clean-dedup', 'id': '661070', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 8346, 'n_tokens_neox': 7438, 'n_words': 2871, 'path': '2019_05_NIH/Day 1/Intro_NumPy.ipynb', 'repo_name': 'mne-tools/mne-workshops', 'size': '21774'}