<jupyter_start><jupyter_text>Lecture 3 Numpy  Covered topics:* `numpy` - a package for storing and working with N-dimensional arrays of data.* Why Python's `list` object is not adequate.  * Dealing with high-dimensional data tricky.  * Does not have 'maths' notation, e.g. $x + y$  * Speed.* Slicing.* Examples:  * Mandelbrot set.  * Markov chains.  * Calculating $\pi$ via random sampling. Working with numbers is at the core of most activities in scientific and engineering computations. The topic is so important that there entire libraries dedicated to working with numbers. `numpy` is the most important and widely used package for dealing with numerical data in Python. We have already seen Python's `list` data structure. We used it for storing a list (sometimes called an array) of data of arbitrary type. We also know how to *slice* the list to get out particular elements we are interested in, and *iterate* over the list.*Reminder:*<jupyter_code>mylist = [0, 3, 4, 5, 'test']
print(mylist[0:3])
print(mylist[-1])
for l in mylist:
    print(l)<jupyter_output>[0, 3, 4]
test
0
3
4
5
test<jupyter_text>'Problems' with Python's `list`.1. Dealing with high-dimension data is tricky.\begin{equation}A = \begin{bmatrix} 0 & 1 & 2 \\3 & 4 & 5 \\6 & 7 & 8\end{bmatrix}\end{equation}Idea: A list of lists! Let's put each row in its own `list` and then put the rows in another list.<jupyter_code>A = [[0, 1, 2], [3, 4, 5], [6, 7, 8]]<jupyter_output><empty_output><jupyter_text>*Exercises*: 1. Select the first row.2. Select the first number of the second row.3. Print the first column. (Not so easy).<jupyter_code>A = [[0, 1, 2], [3, 4, 5], [6, 7, 8]]

print(A[0])
print(A[1][0])
for row in A:
    print(row[0])<jupyter_output>[0, 1, 2]
3
0
3
6<jupyter_text>2. `list` does not conform to 'standard' mathematical conventions.$$ x = [1, 2, 3]$$$$ y = [4, 5, 6]$$$$ x + y = \; ? $$<jupyter_code>x = [1, 2, 3]
y = [4, 5, 6] 
x + y<jupyter_output><empty_output><jupyter_text>Something better: `numpy`<jupyter_code>import numpy
import numpy as np # np is now an alias (shorthand) for numpy
np.lookfor('create array')
x = np.array([1, 2, 3])
y = np.array([4, 5, 6])
# f-strings, Python 3.6 and above
print(f"x: {x}")
print(f"y: {y}")
x + y<jupyter_output><empty_output><jupyter_text>3. Speed *Exercise*: Write a function `numpy_add` that takes two numpy arrays, loops over the entries, and computes the sum entry-by-entry and places the result into a third numpy array.Tips:* `def` - define a function* `return x` - return `x` from the function* `len(a)` - get the length of an array.* `np.zeros(n)` - make a numpy array of length `n`* `for i in range(0, len(a)): print(i)` - counting indices $0, ..., \mathrm{len}(a) - 1$<jupyter_code>def numpy_add(a, b):
    assert(len(a) == len(b))
    c = np.zeros(len(a))
    for i in range(0, len(a)):
        c[i] = a[i] + b[i]
    return c
        
a = np.ones(1000000)
b = np.ones(1000000)
%time c = numpy_add(a, b)
%time c = a + b<jupyter_output>CPU times: user 8.84 ms, sys: 2.78 ms, total: 11.6 ms
Wall time: 10.8 ms<jupyter_text>Key points:* Avoid writing 'tight' loops in Python.  * Tight loop - looping over thousands of entries of data and performing a single operation on each piece of data.* Think about how algorithms can be written without tight loops.* Always prefer functionality in numpy that can act on many pieces of data in one call. *Example*: Plot the function $\sin(x)$.<jupyter_code>x = np.linspace(0.0, 2.0*np.pi, num=50)
y = np.sin(x) # Apply sin function to all of the entries of x
print(f"x: {x}")
print(f"y: {y}")
%matplotlib inline
import matplotlib.pyplot as plt
plt.xlabel(r"$x$")
plt.ylabel(r"$\sin(x)$")
plt.plot(x, y, 'o-')<jupyter_output><empty_output><jupyter_text>Slicing *Slicing* allows us to select one or more than one elements of an existing array.<jupyter_code>a = np.arange(5)
print(a)<jupyter_output>[0 1 2 3 4]<jupyter_text>Selecting single elements of an array * Select the first element<jupyter_code>a[0]<jupyter_output><empty_output><jupyter_text>* Select the third element (remember, index starts from zero!)<jupyter_code>a[2]<jupyter_output><empty_output><jupyter_text>* Select the fifth element (remember, index starts from zero!)<jupyter_code>a[4]<jupyter_output><empty_output><jupyter_text>* Select the last element, without knowing the length of the array.<jupyter_code>a[-1]<jupyter_output><empty_output><jupyter_text>Selecting multiple elements of an array  * Select the first and second elements of the array.* syntax: `[start:end]` exclusive of `end`  * default start: `0`  * default end: `N`<jupyter_code>print(a[0:len(a)])
print(a[:])<jupyter_output>[0 1 2 3 4]
[0 1 2 3 4]<jupyter_text>* First to end<jupyter_code>a[0:5]
a[0:]<jupyter_output><empty_output><jupyter_text>* Exercise: Select the second, third and fourth elements<jupyter_code>a[1:4]<jupyter_output><empty_output><jupyter_text>* Exercise: Select everything.<jupyter_code>a[:]<jupyter_output><empty_output><jupyter_text>Advanced slicing * Select every second element.* Syntax: `[start:end:stride]`  * default start: `0` for `stride > 0` and `-1` for `stride < 0`  * default end: `N` for `stride > 0` and `- N - 1` for `stride < 0`   * default stride: `1`  Full details: https://numpy.org/doc/stable/reference/arrays.indexing.html<jupyter_code>a[::2]
a[::2]<jupyter_output><empty_output><jupyter_text>Let's make a longer array.<jupyter_code>b = np.arange(0, 10)
print(b)<jupyter_output>[0 1 2 3 4 5 6 7 8 9]<jupyter_text>* Exercise: Select every third entry starting with the second entry.<jupyter_code>b[1::3]<jupyter_output><empty_output><jupyter_text>* *Exercise*: Starting with the last entry, select every third entry working backwards through the array. (Hint: Try playing with a `stride` of `-1`).<jupyter_code>b[-1::-3]<jupyter_output><empty_output><jupyter_text>* *Exercise*: Starting with the second to last entry, select every third entry working backwards through the array.<jupyter_code>b[-2::-3]<jupyter_output><empty_output><jupyter_text>N-dimensional indexing* Slices in each axis are seperated by commas.* `[start0:stop0:step0, start1:stop1:step1]`  * Default: `[,]` = `[::,::]`<jupyter_code>c = np.arange(0, 32).reshape(8, 4)
print(c)<jupyter_output>[[ 0  1  2  3]
 [ 4  5  6  7]
 [ 8  9 10 11]
 [12 13 14 15]
 [16 17 18 19]
 [20 21 22 23]
 [24 25 26 27]
 [28 29 30 31]]<jupyter_text>* Select the element on the third row in the fourth column.<jupyter_code>c[2, 3] # [axis0, axis1]<jupyter_output><empty_output><jupyter_text>* *Exercise*: From every second row, select every element.  * Select every second row.<jupyter_code>c[::2, ::]
c[::2, ]
c[::2]<jupyter_output><empty_output><jupyter_text>* Exercise: Select the elements in every third row, starting with the second row, and in every second column.<jupyter_code>c[1::3, ::2]<jupyter_output><empty_output><jupyter_text>Slicing based on integers* Select the elements at (1, 1), (2, 2) and (3, 3).<jupyter_code>c[(1, 2, 3), (1, 2, 3)]<jupyter_output><empty_output><jupyter_text>* Exercise: Select the elements at (5, 3) and (1, 2).<jupyter_code>c[[5, 1], [3, 2]] # Can use lists, ndarrays or tuples.<jupyter_output><empty_output><jupyter_text>Slicing based on booleans<jupyter_code>c < 5
c < 5<jupyter_output><empty_output><jupyter_text>Mini-project: Mandelbrot set Key Skills:* Complex numbers* Creating a grid of points in the complex plane* Iteration* Avoiding tight loops with operations on `numpy` arrays* Boolean masking* Basic plot The Mandelbrot set is the set of complex numbers $c$ that remain bounded under repeated application of the function:$$z_0 = c \quad z_{n + 1} = z_{n}^2 + c$$  A reminder of complex numbers$$z = x + j y$$ where $j = \sqrt{-1}$ and $x = \Re(z)$ and $y = \Im(z)$ are the real and complex parts of $z$.<jupyter_code>z = 3.0 + 1j*4.0
print("Real part of z: {}".format(np.real(z)))
print("Complex part of z: {}".format(np.complex(z)))
print("z + z = {}".format(z + z))
print("z**2 = {}".format(z**2))<jupyter_output>Real part of z: 3.0
Complex part of z: (3+4j)
z + z = (6+8j)
z**2 = (-7+24j)<jupyter_text>Step 1: Create a grid of points $c$ in the complex plane1. Tool: `np.linspace`. Create 300 points between -2 and 1 and place the result in `x`.2. Tool: `np.linspace`. Create 300 points between -1.5 and 1.5 and place the result in `y`.3. Tool: `x[:, np.newaxis]`, `y[np.newaxis, :]`. What does this do?<jupyter_code>x = np.linspace(-2, 1, 300)
y = np.linspace(-1.5, 1.5, 300)
c = x[:, np.newaxis] + 1j*y[np.newaxis, :]<jupyter_output><empty_output><jupyter_text>Step 2: Apply the function The Mandelbrot set is the set of complex numbers $c$ that remain bounded under repeated application of the function:$$z_0 = c \quad z_{n + 1} = z_{n}^2 + c$$1. Tool: `for i in range(0, 50):`2. `numpy` element-wise operations `z**2` `+` `c`3. Making a copy of c `np.zeros_like(c)`<jupyter_code>num_iterations = 50
z = np.zeros_like(c)
for _ in range(0, num_iterations):
    z = z**2 + c<jupyter_output><empty_output><jupyter_text>Step 3: Select the $c$ "that remain bounded"1. Tool: Boolean operation `z < bound`<jupyter_code>bound = 100.0
mandelbrot_set = z < bound<jupyter_output><empty_output><jupyter_text>Step 4: Plot the points $c$ that are in the Mandelbrot set1. Tool: `plt.imshow`.<jupyter_code>import matplotlib.pyplot as plt
plt.imshow(mandelbrot_set.T, extent=[-2, 1, -1.5, 1.5])
plt.gray()
plt.xlabel("$\Re(c)$")
plt.ylabel("$\Im(c)$")
plt.show()<jupyter_output><empty_output><jupyter_text>Mini-project: Markov ChainKey Skills:* Calculating a transition matrix* Defining a transition matrix in `numpy`* Matrix-vector multiplication Wikipedia: A *Markov chain* is a *stochastic model* describing a sequence of possible events in which the *probability of each event* depends only on the state attained in the previous event. Imagine you are standing in a five room house.*Drawing of floorplan of house*The doors between the rooms of the house are locked.Problem: Every hour the doors of the house are unlocked, and *you must* move to another room in the house. You have no preference over which room you want to be in next, so you move randomly between the rooms. **Which room will you spend the most time in?** Now imagine you are in room 1. You could move to rooms 2, 4 and 5 (three rooms). You cannot move to room 3.Denote $m_{ij}$ the *probability* of moving between room $j$ and room $i$.What are $m_{21}$, $m_{41}$ and $m_{51}$? What are $m_{11}$ and $m_{31}$? $$m_{21} = m_{41} = m_{51} = 1/3$$$$m_{11} = m_{31} = 0$$ ${m}_{1} = \begin{bmatrix}  m_{11} \\ m_{21} \\ m_{31} \\ m_{41}  \\m_{51} \end{bmatrix} = \begin{bmatrix}  0 \\ 1/3 \\ 0 \\ 1/3  \\ 1/3 \end{bmatrix} $ *Exercise*: Write down $m_1$, $m_2$, $m_3$, $m_4$, $m_5$. The transition matrix $M$ has entries $m_{ij}$.*Exercise*: Write the transition matrix $M$ on paper.*Exercise*: Write the transition matrix $M$ in a `numpy` array. $$M = \begin{bmatrix} 0 & 1/3 & 0 & 1/3 & 1/4 \\1/3 & 0 & 1/3 & 0 & 1/4 \\0 & 1/3 & 0 & 1/3 & 1/4 \\1/3 & 0 & 1/3 & 0 & 1/4 \\1/3 & 1/3 & 1/3 & 1/3 & 0 \\\end{bmatrix}$$<jupyter_code>M = np.array([[0, 1/3, 0, 1/3, 1/4], 
              [1/3, 0, 1/3, 0, 1/4], 
              [0, 1/3, 0, 1/3, 1/4],
              [1/3, 0, 1/3, 0, 1/4],
              [1/3, 1/3, 1/3, 1/3, 0]])
print(M)<jupyter_output>[[0.         0.33333333 0.         0.33333333 0.25      ]
 [0.33333333 0.         0.33333333 0.         0.25      ]
 [0.         0.33333333 0.         0.33333333 0.25      ]
 [0.33333333 0.         0.33333333 0.         0.25      ]
 [0.33333333 0.33333333 0.33333333 0.33333333 0.        ]]<jupyter_text>Question: Which room will you spend the most time in?Let's assume we are certain we start in room 1.$$ p_0 = \begin{bmatrix} 1.0 & 0.0 & 0.0 & 0.0 & 0.0 \end{bmatrix} $$<jupyter_code>p_0 = np.array([1.0, 0.0, 0.0, 0.0, 0.0])<jupyter_output><empty_output><jupyter_text>Where will we end up next?$$ p_1 = M p_0$$<jupyter_code>p_1 = M@p_0
print(p_1)
p_2 = M@p_1
print(p_2)
p_3 = M@p_2
print(p_3)
p = np.array([1.0, 0.0, 0.0, 0.0, 0.0])
for i in range(0, 50):
    p = M@p
print(p)<jupyter_output>[0.1875 0.1875 0.1875 0.1875 0.25  ]<jupyter_text>*Question*: What happens if we are not sure where we start?<jupyter_code>p = np.array([0.25, 0.25, 0.0, 0.25, 0.25])
for i in range(0, 50):
    p = M@p
print(p)<jupyter_output>[0.1875 0.1875 0.1875 0.1875 0.25  ]<jupyter_text>Mini-project: Calculating $\pi$ with Monte CarloKey skills:* Simulating random numbers* Boolean operations* Slicing with strides* Counting* Algorithms with no tight-loops Imagine you are blind-folded and you have an infinite number of darts to throw at the following square board:*Draw square with quarter circle*How can we calculate the value of $\pi$?*Monte Carlo* (random simulation) provides a powerful answer to this problem. Steps:    1. Randomly throw (many) darts into the square (random simulation).2. Calculate the ratio $r$ of the number of darts that fall inside the circle/total number thrown.3. $\pi = 4r$  Step 1: Simulate from uniform distributionWe want to 'simulate' the idea of throwing darts randomly at a square board.Computers can *generate* (pseudo-)random numbers.1. Tool: `np.random.uniform(size=())` to simulate N draws of uniformly distributed random numbers on the board.Uniformly distributed - equally likely to be at any value between $[0, 1]$.<jupyter_code>N = 100000
darts = np.random.uniform(size=(N, 2))
plt.axis('equal')
plt.xlabel('x')
plt.ylabel('y')
plt.plot(darts[::100, 0], darts[::100, 1], 'x')<jupyter_output><empty_output><jupyter_text>Step 2: Calculate the squared distance from the point to the origin.1. Tool: Slicing to get all x and all y coordinates. 2. Tool: `**2` to square entries of an array.<jupyter_code>dist = darts[:, 0]**2 + darts[:, 1]**2<jupyter_output><empty_output><jupyter_text>Step 3: Which darts landed inside the circle?1. Tool: Boolean selection.<jupyter_code>inside_outside = dist < 1.0
print(inside_outside)<jupyter_output>[ True  True  True ...  True False  True]<jupyter_text>Step 4: Count the number of `True` values.1. Tool: Use `np.count_nonzero`.<jupyter_code>np.lookfor('count')<jupyter_output>Search results for 'count'
--------------------------
numpy.bincount
    Count number of occurrences of each value in array of non-negative ints.
numpy.busday_count
    Counts the number of valid days between `begindates` and
numpy.nanmax
    Return the maximum of an array or maximum along an axis, ignoring any
numpy.nanmin
    Return minimum of an array or minimum along an axis, ignoring any NaNs.
numpy.count_nonzero
    Counts the number of non-zero values in the array ``a``.
numpy.ma.count
    Count the non-masked elements of the array along the given axis.
numpy.fromfile
    Construct an array from data in a text or binary file.
numpy.fromiter
    Create a new 1-dimensional array from an iterable object.
numpy.nancumsum
    Return the cumulative sum of array elements over a given axis treating Not a
numpy.bytes0.count
    Return the number of non-overlapping occurrences of subsection sub in
numpy.frombuffer
    Interpret a buffer as a 1-dimensional array.
numpy.fromstring
    A new[...]<jupyter_text>2. Found tool: np.count_nonzero<jupyter_code>num_inside = np.count_nonzero(inside_outside)
print(num_inside)<jupyter_output>78599<jupyter_text>Step 5. Calculate $\pi$!<jupyter_code>estimate_pi = 4.0*num_inside/N
print(estimate_pi)
print(np.pi)
print(np.abs(np.pi - estimate_pi))<jupyter_output>3.14396
3.141592653589793
0.00236734641020675
{'fraction_non_alphanumeric': 0.1421457673155274, 'fraction_numerical': 0.05238310806441182, 'mean_word_length': 4.475920679886686, 'pattern_counts': {'":': 0, '<': 165, '<?xml version=': 0, '>': 160, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 12, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '43', 'dir': 'jupyter-structured-clean-dedup', 'id': '267351', 'license': 'MIT', 'n_tokens_mistral': 6165, 'n_tokens_neox': 5616, 'n_words': 1879, 'path': 'notebooks/03_numpy.ipynb', 'repo_name': 'jhale/scientific-python', 'size': '15463'}