<jupyter_start><jupyter_text>NumPy -  Multidimensional Arrays The numpy module is an extremely powerful package for dealing with numerical calculations.  At its core it provides a way for users to create and manipulate N-dimensional arrays.  It also provides a number of mathematical functions (e.g. exponential and trigonometric functions) and solvers (e.g. FFTW).More importantly, these arrays behave much like arrays in C, C++, and Fortran:* Arrays have a fixed size at creation (no dynamic growth)* Homogeneous datatypes (e.g. all floats).  Means known, fixed memory size for array* Contiguous in memory (improved performance)These features allow for better performance, while still maintaining the flexibility and ease of use of Python. To get started, we'll import the NumPy module:<jupyter_code>import numpy as np<jupyter_output><empty_output><jupyter_text>Creating NumPy ArraysYou can initialise NumPy arrays in a variety of ways:* Python lists and tuples* NumPy functions (`arange`, `zeros`, etc.)* Input data from files  Lists/tuples<jupyter_code># vector example
x = np.array([1,2,3,4])
x
A = np.array([[1, 2], [3, 4]])
A<jupyter_output><empty_output><jupyter_text>We can inspect the type of the arrays:<jupyter_code>type(x), type(A)<jupyter_output><empty_output><jupyter_text>If we want to know more about the dimensions of our arrays, we can use the `size` and `shape` functions:<jupyter_code>x.shape
A.shape
A.size<jupyter_output><empty_output><jupyter_text>Note that these functions are acutally called as `ndarray.size` and `ndarray.shape`.  They are functions available to ndarray objects via the NumPy module.  Equivalently, you can call the same functions directly from the NumPy module:<jupyter_code>np.shape(A)<jupyter_output><empty_output><jupyter_text>These arrays look very much like Python lists...so why not just use a list for computation?  Python is great with lists, why do we need a new datatype?Performance.As mentioned in the slides, Python lists are dynamically typed, have an inefficient memory layout (at least for numerical computations), and don't support common array operations like dot products.NumPy arrays only support certain datatypes:<jupyter_code>A[0,0] = "Brian"
A.dtype<jupyter_output><empty_output><jupyter_text>Supported datatypes include `int`, `float`, `complex`, `bool`, Python objects, and more.  We can even set bit size for these:<jupyter_code>B = np.array([[1,3],[2,4]], dtype=complex)
B<jupyter_output><empty_output><jupyter_text>A note on memory layoutWe've mentioned that NumPy stores arrays in contiguous memory.  But what does that mean, and how does it apply to multi-dimensional arrays?  Simply, a multidimensional array is actually "stretched out" and stored in memory as a long, 1D array.  There are a few ways we can "stretch" our arrays, and they're known as ***row major*** and ***column major*** order:  Why does this matter?Imagine we traverse this 2D array in a for loop, and assume it's stored in row major form, where `n` is the number of rows and `m` is the number of columns.  What happens if we access the code like this? `for j in 0:m    for i in 0:n        A[i][j] += x[i]` The problem is we'll be jumping around in memory, destroying our performance.  Modern compilers try to pull in contiguous blocks of data from RAM, and put them in cache (smaller, but faster, regions of memory closer to the CPU).  This staging process provides a LOT of performance, and choosing the wrong storage format (or using the wrong indices) will destroy your performance. By default, NumPy stores arrays in row major order:<jupyter_code>A = np.random.rand(3,3)
A.flags<jupyter_output><empty_output><jupyter_text>The hint here is ***C_CONTIGUOUS***.  Traditionally, C/C++ store data in row major form, and Fortran stores it in column major form.  The convention carries over still today, and so Python denotes row major as `C` and column major as `F` Python actually allows you to change the order (if you want):<jupyter_code>A = np.random.rand(3,3)
A = np.asfortranarray(A)
A.flags<jupyter_output><empty_output><jupyter_text>Functions Generally we don't manually set up the elements of an array; instead we can use built-in functions to generate arrays.  arange<jupyter_code># start, stop, step size
x = np.arange(0, 5, 1)  
x

x = np.arange(-1, 1, 0.2)  
x<jupyter_output><empty_output><jupyter_text>mgridSimilar to `meshgrid` in MATLAB, it generates 2 arrays where the values correspond to the indices.<jupyter_code>x, y = np.mgrid[0:4, 0:4]
x
y<jupyter_output><empty_output><jupyter_text>rand NumPy also allows you to generate random numbers:<jupyter_code># Generates uniform range of numbers in [0,1]
A = np.random.rand(3,3)
A<jupyter_output><empty_output><jupyter_text>zeros and ones<jupyter_code>x = np.zeros(3)
x
x = np.ones((5,5))
x<jupyter_output><empty_output><jupyter_text>I/OWe can read in CSV data into a NumPy array as well.  We'll use some example temperature data along with the `genfromtxt` function:<jupyter_code>!head ../example_data/temp_data.dat
import numpy as np
data = np.genfromtxt('../demos/data/temp_data.dat')
data.shape
import matplotlib.pyplot as plt
%matplotlib inline
fig, ax = plt.subplots(figsize=(14,4))
ax.plot(data[:,0]+data[:,1]/12.0+data[:,2]/365, data[:,5])
ax.axis('tight')
ax.set_title('Temperatures in Stockholm')
ax.set_xlabel('Year')
ax.set_ylabel('Temperature (C)')<jupyter_output><empty_output><jupyter_text>We can also write out CSV files from NumPy arrays:<jupyter_code>A = np.random.rand(5,5)
A
np.savetxt("../demos/data/rand_mat.csv",A)
!head ../demos/data/rand_mat.csv
# Let's specify a file format
np.savetxt("../demos/data/rand_mat.csv", A, fmt='%.3f')
!cat ../demos/data/rand_mat.csv<jupyter_output>0.410 0.652 0.450 0.565 0.477
0.767 0.007 0.833 0.163 0.809
0.152 0.424 0.661 0.158 0.406
0.849 0.727 0.263 0.081 0.935
0.917 0.611 0.202 0.068 0.173<jupyter_text>Manipulating arrays: indexing Indexing arrays in NumPy is done via square brackets:<jupyter_code>A[1,1]
A[1][1]<jupyter_output><empty_output><jupyter_text>With multi-dimensional arrays, if you omit an index it will show the whole row:<jupyter_code>A[1]<jupyter_output><empty_output><jupyter_text>Assigning values works the same:<jupyter_code>A[0,0] = 0.5
A<jupyter_output><empty_output><jupyter_text>We can also use `:` to access all elements in a row or column (array slicing):<jupyter_code>A[1,:]
A[:,1]
A[0,:] = 1
A
# Slice out several columns of A
cols = A[:,[0,2]]
cols<jupyter_output><empty_output><jupyter_text>***Note:*** These slices aren't copies...they just point to the original array<jupyter_code># A before we modify it
A
# Pull out row 3 from A
test_row = A[3,:]
test_row
# Zero out the row
test_row[:] = 0
test_row
# Print A
A<jupyter_output><empty_output><jupyter_text>If we want an explicit copy, we can use the `copy` function:<jupyter_code>test_row = A[2,:]
copyrow = test_row.copy()
copyrow<jupyter_output><empty_output><jupyter_text>Linear algebra NumPy allows for us for to do common operations like matrix-vector multiplication, scalar multiplications, dot products, etc.  These operations are so important to computation that entire libraries have been developed and optimized just to handle things like `a*x+y` or `A*B=C` (e.g. BLAS, LAPACK, MKL, etc.).  NumPy is able to provide performant versions of these operations via something called ***ufuncs*** (universal functions).  Ufuncs provide a convenient interface to ***vectorized*** routines in compiled libraries.<jupyter_code>def compute_square(vals):
    sq_vals = np.zeros(len(vals))
    for i in range(len(vals)):
        sq_vals[i] = vals[i] * vals[i]
    return sq_vals

compute_square(np.random.randint(1, 10, size=5))
big_array = np.random.randint(1, 100, size=1000000)
%timeit compute_square(big_array)<jupyter_output>652 ms ± 47 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)<jupyter_text>Imagine you have a 2D grid of points, ~3000 points in each direction....if you have to compute a square for each value (which is not an uncommon operation) you'll be waiting several seconds to do this.  Now imagine you have to this several thousand (or million) times in a simulation.  As we've mentioned earlier, the delay is not the multiplication operation itself, but rather all of the type checking and dynamic lookups the Python interpreter has to do.Now let's try a ufunc:<jupyter_code>%timeit np.square(big_array)<jupyter_output>844 µs ± 7.79 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)<jupyter_text>Massive difference.  Element-wise Operations<jupyter_code>x = np.arange(10)
x
# Scalar addition
x+1
# Exponentiation
2**x
# Division
a = np.random.randint(1,100,size=10000)
b = np.random.randint(1,100,size=10000)

%timeit c = a / b<jupyter_output>36.1 µs ± 1 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)<jupyter_text>Note that each of these operations is simply a wrapper around a specific NumPy function: | Operation | NumPy Function | Action || --- | --- | --- || `+` | `np.add` | Addition (e.g., 1 + 1 = 2) ||`-`|`np.subtract`|Subtraction (e.g., 3 - 2 = 1)||`-`|`np.negative`|Unary negation (e.g., -2)||`*`|`np.multiply`|Multiplication (e.g., 2 * 3 = 6)||`/`|`np.divide`|Division (e.g., 3 / 2 = 1.5)||`//`|`np.floor_divide`|Floor division (e.g., 3 // 2 = 1)||`**`|`np.power`|Exponentiation (e.g., 2 ** 3 = 8)||`%`|`np.mod`|Modulus/remainder (e.g., 9 % 4 = 1)|  Dot productBe careful if you're trying to do things like dot products.  Remember, the above functions are **element-wise** operations:<jupyter_code>a.shape
b.shape
a * b<jupyter_output><empty_output><jupyter_text>Normally, we'd expect a scalar value (remember, a dot product between vectors produces a scalar value).  We can do this with the `dot` function:<jupyter_code>np.dot(a, b)
# Example with matrices
A = np.random.randint(1,10,size=(5,5))
B = np.random.randint(1,10,size=(5,5))

# Matrix-vector product
x = np.random.rand(5)
np.dot(A,x)
np.dot(A,B)<jupyter_output><empty_output><jupyter_text>We can also cast these arrays as matrices, using the `matrix` function.  This lets us use the arithmetic operators for **matrix algebra**:<jupyter_code># Cast a new matrix
M = np.matrix(A)

# Create a column vector
# Note: T is the transpose operator
v = np.matrix(x).T

M * v<jupyter_output><empty_output><jupyter_text>Matrix ComputationsNumPy's `linalg` module provides a range of common operations for computing things like inverses, norms, eigenvalues, decompositions, and even basic linear solvers. Inverse<jupyter_code>np.linalg.inv(A)
# Equivalent matrix function
C = np.matrix(A)

C.I<jupyter_output><empty_output><jupyter_text>Just to check:<jupyter_code>I_mat = C.I * C
np.set_printoptions(precision=1)

print(I_mat)<jupyter_output>[[ 1.0e+00 -2.2e-16 -8.9e-16  2.2e-16  4.4e-16]
 [ 1.9e-16  1.0e+00 -5.6e-17  4.2e-17  1.9e-16]
 [ 1.1e-16  0.0e+00  1.0e+00  2.8e-16  1.1e-16]
 [ 4.4e-16  0.0e+00 -1.3e-15  1.0e+00  4.4e-16]
 [ 4.4e-16  0.0e+00  8.9e-16  2.2e-16  1.0e+00]]<jupyter_text>Determinants<jupyter_code>np.linalg.det(C)<jupyter_output><empty_output><jupyter_text>Eigenvalues<jupyter_code>w,v = np.linalg.eig(A)
w
v<jupyter_output><empty_output><jupyter_text>Returns an object with 2 arrays: the first array are the eigenvalues, and the second array is composed of the eigenvectors (normalised).Note that the `ith` column of v (`v[:,i]`) corresponds to the `ith` eigenvalue of w (`w[i]`)  NormsNumPy provides a lot of different norms to choose from:<jupyter_code># 2-norm for vectors
np.linalg.norm(x)
# 2-norm for matrices (aka Frobenius)
np.linalg.norm(A)
# Infinity norm for vectors
np.linalg.norm(x,np.inf)
# Negative infinity norm for matrices
np.linalg.norm(A, -np.inf)<jupyter_output><empty_output><jupyter_text>Max and Min<jupyter_code>A.min()
A
x.max()<jupyter_output><empty_output><jupyter_text>Reshaping Arrays NumPy has methods for quickly and efficiently manipulating arrays that don't involve making a copy of the data (which would greatly impact performance)<jupyter_code>A
n, m = A.shape

# Flatten the matrix into an n*m 1D array
B = A.reshape((1,n*m))
B<jupyter_output><empty_output><jupyter_text>Let's alter the array:<jupyter_code>B[0,0:10] = 0
B
A<jupyter_output><empty_output><jupyter_text>Original values in A were changed...B is simply a pointer to AThere are situations where you want to reshape and make a copy.  In that case, you can use the `flatten` function:<jupyter_code>B = A.flatten()
B[0,0:10] = -1
B[0:10] = -1
B
A<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.11624548736462094, 'fraction_numerical': 0.032089851584436425, 'mean_word_length': 5.34725050916497, 'pattern_counts': {'":': 0, '<': 140, '<?xml version=': 0, '>': 140, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '36', 'dir': 'jupyter-structured-clean-dedup', 'id': '148956', 'license': 'CC0-1.0', 'n_tokens_mistral': 4629, 'n_tokens_neox': 4334, 'n_words': 1435, 'path': 'notebooks/01-NumPy-Intro.ipynb', 'repo_name': 'PawseySC/Using-Python-in-HPC', 'size': '12465'}