<jupyter_start><jupyter_text>Why do we need numpy?* You may have heard "Python is slow", this is true when it concerns looping over many small python objects* Python is dynamically typed and everything is an object, even an `int`. There are no primitive types.* Numpy's main feature is the `ndarray` class, a fixed length, homogeniously typed array class.* Numpy implements a lot of functionality in fast c, cython and fortran code to work on these arrays* python with vectorized operations using numpy can be blazingly fastSee: [Python is not C](https://www.ibm.com/developerworks/community/blogs/jfp/entry/Python_Is_Not_C?lang=en)But the most important reason:* More beautiful code<jupyter_code>import numpy as np<jupyter_output><empty_output><jupyter_text>More beautiful code through vectorisationpure python with list comprehension<jupyter_code>voltages = [10.1, 15.1, 9.5]
currents = [1.2, 2.4, 5.2]

resistances = [U * I for U, I in zip(voltages, currents)]
resistances<jupyter_output><empty_output><jupyter_text>Using numpy<jupyter_code>U = np.array([10.1, 15.1, 9.5])
I = np.array([1.2, 2.4, 5.2])

R = U * I
R<jupyter_output><empty_output><jupyter_text>Finding the point with the smallest distance<jupyter_code>import math

def euclidean_distance(p1, p2):
    return math.sqrt((p1[0] - p2[0])**2 + (p1[1] - p2[1])**2)

point = (1, 2)
points = [(3, 2), (4, 2), (3, 0)]

min_distance = float('inf')
for other in points:
    distance = euclidean_distance(point, other)
    if distance < min_distance:
        closest = other
        min_distance = distance 

print(min_distance, closest)
point = np.array([1, 2])
points = np.array([(3, 2), (4, 2), (3, 0)])

distance = np.linalg.norm(point - points, axis=1)
idx = np.argmin(distance)

print(distance[idx], points[idx])<jupyter_output>2.0 [3 2]<jupyter_text>Small example timings<jupyter_code>import math


def var(data):
    '''
    knuth's algorithm for one-pass calculation of the variance
    Avoids rounding errors of large numbers when doing the naive
    approach of `sum(v**2 for v in data) - sum(v)**2`
    '''
    
    n = 0
    mean = 0.0
    m2 = 0.0
    
    if len(data) < 2:
        return float('nan')

    for value in data:
        n += 1
        delta = value - mean
        mean += delta / n
        delta2 = value - mean
        m2 += delta * delta2

    return m2 / n 
list(range(10))
%%timeit

l = list(range(1000))
var(l)
%%timeit

a = np.arange(1000)  # array with numbers 0,...,999

np.var(a)<jupyter_output>31.1 µs ± 7.2 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)<jupyter_text>Basic math: vectorizedOperations on numpy arrays work vectorized, element-by-element** Lose your loops **<jupyter_code># create a numpy array from a python a python list
a = np.array([1.0, 3.5, 7.1, 4, 6])
2 * a
a**2
a**a
np.cos(a)<jupyter_output><empty_output><jupyter_text>**Attention: You need the `cos` from numpy!**<jupyter_code>math.cos(a)<jupyter_output><empty_output><jupyter_text>Most normal python functions with basic operators like `*`, `+`, `**` simply work becauseof operator overloading:<jupyter_code>def poly(x):
    return x + 2 * x**2 - x**3

poly(a)
poly(np.e), poly(np.pi)<jupyter_output><empty_output><jupyter_text>Useful properties<jupyter_code>len(a)
a.shape
a.dtype
a.ndim<jupyter_output><empty_output><jupyter_text>Arbitrary dimension arrays<jupyter_code># two-dimensional array
y = np.array([[1, 2, 3],
              [4, 5, 6],
              [7, 8, 9]])

y + y
## since python 3.5 @ is matrix product
y @ y
# Broadcasting, changing array dimensions to fit the larger one

y + np.array([1, 2, 3])<jupyter_output><empty_output><jupyter_text>Reduction operationsNumpy has many operations, which reduce dimensionality of arrays<jupyter_code>x = np.random.normal(0, 1, 1000)
np.sum(x)
np.prod(x)
np.mean(x)<jupyter_output><empty_output><jupyter_text>Standard Deviation<jupyter_code>np.std(x)<jupyter_output><empty_output><jupyter_text>Standard error of the mean<jupyter_code>np.std(x, ddof=1) / np.sqrt(len(x))<jupyter_output><empty_output><jupyter_text>Sample Standard Deviation<jupyter_code>np.std(x, ddof=1)<jupyter_output><empty_output><jupyter_text>Most of the numpy functions are also methods of the array<jupyter_code>x.mean(), x.std(), x.max(), x.min()<jupyter_output><empty_output><jupyter_text>Difference between neighbor elements<jupyter_code>z = np.arange(10)**2
diff_z = np.diff(z)

print(z)
print(diff_z)<jupyter_output>[ 0  1  4  9 16 25 36 49 64 81]
[ 1  3  5  7  9 11 13 15 17]<jupyter_text>Reductions on multi-dimensional arrays<jupyter_code>array2d = np.arange(20).reshape(4, 5)

array2d
np.sum(array2d, axis=0)
np.var(array2d, axis=1)<jupyter_output><empty_output><jupyter_text>Exercise 1Write a function that calculates the analytical linear regression for a set ofx and y values.Reminder:$$ f(x) = a \cdot x + b$$with $$\hat{a} = \frac{\mathrm{Cov}(x, y)}{\mathrm{Var}(x)} \\\hat{b} = \bar{y} - \hat{a} \cdot \bar{x}$$<jupyter_code># %load 04_01_numpy_solutions/exercise_linear.py
x = np.linspace(0, 1, 50)
y = 5 * np.random.normal(x, 0.1) + 2  # see section on random numbers later

a, b = linear_regression(x, y)
a, b<jupyter_output><empty_output><jupyter_text>Helpers for creating arrays<jupyter_code>np.zeros(10)
np.ones((5, 2))
np.full(5, np.nan)
np.empty(10)  # attention, uninitialised memory, be carefull
np.linspace(-2, 1, 1)
# like range() for arrays:
np.arange(5)
np.arange(2, 10, 2)
np.logspace(-4, 5, 10)
np.logspace(1, 4, 4, base=2)<jupyter_output><empty_output><jupyter_text>Numpy Indexing* Element access* Slicing<jupyter_code>x = np.arange(0, 10)

# like lists:
x
# like lists:
x[0]
# all elements with indices ≥1 and <4:
x[1:4]
# negative indices count from the end
x[-1], x[-2]
# combination:
x[3:-2]
# step size
x[::2]
# trick for reversal: negative step
x[::-1]
y = np.array([x, x + 10, x + 20, x + 30])
y
# only one index ⇒ one-dimensional array
y[2]
# other axis: (: alone means the whole axis)
y[:, 3]
# inspecting the number of elements per axis:
y[:, 1:3].shape<jupyter_output><empty_output><jupyter_text>Changing array content<jupyter_code>y
y[:, 3] = 0
y<jupyter_output><empty_output><jupyter_text>Using slices on both sides<jupyter_code>y[:,0] = x[3:7]
y<jupyter_output><empty_output><jupyter_text>Transposing inverts the order of the dimensions<jupyter_code>y
y.shape
y.T
y.T.shape<jupyter_output><empty_output><jupyter_text>Masks* A boolean array can be used to select only the element where it contains `True`.* Very powerfull tool to select certain elements that fullfill a certain condition<jupyter_code>a = np.linspace(0, 2, 11)
b = np.random.normal(0, 1, 11)

print(b >= 0)
print(a[b >= 0])
a[[0, 2]] = np.nan
a
a[np.isnan(a)] = -1
a<jupyter_output><empty_output><jupyter_text>Random numbers* numpy has a larger number of distributions builtin<jupyter_code>np.random.uniform(-1, 1, 10)
np.random.normal(0, 5, (2, 10))
np.random.poisson(5, 2)<jupyter_output><empty_output><jupyter_text>Calculating pi through monte-carlo simulation* We draw random numbers in a square with length of the sides of 2* We count the points which are inside the circle of radius 1The area of the square is$$A_\mathrm{square} = a^2 = 4$$The area of the circle is$$A_\mathrm{circle} = \pi r^2 = \pi$$With $$\frac{n_\mathrm{circle}}{n_\mathrm{square}} = \frac{A_\mathrm{circle}}{A_\mathrm{square}}$$We can calculate pi:$$\pi = 4 \frac{n_\mathrm{circle}}{n_\mathrm{square}}$$<jupyter_code>n_square = 10000000

x = np.random.uniform(-1, 1, n_square)
y = np.random.uniform(-1, 1, n_square)

radius = np.sqrt(x**2 + y**2)

n_circle = np.sum(radius <= 1.0)

print(4 * n_circle / n_square)<jupyter_output>3.1424804<jupyter_text>Exercise1. Draw 10000 gaussian random numbers with mean of $\mu = 2$ and standard deviation of $\sigma = 3$2. Calculate the mean and the standard deviation of the sample3. What percentage of the numbers are outside of $[\mu - \sigma, \mu + \sigma]$?4. How many of the numbers are $> 0$?5. Calculate the mean and the standard deviation of all numbers ${} > 0$<jupyter_code># %load 04_01_numpy_solutions/exercise_gaussian.py<jupyter_output><empty_output><jupyter_text>ExerciseMonte-Carlo uncertainty propagation* The hubble constant as measured by PLANCK is$$H_0 = (67.74 \pm 0.47)\,\frac{\mathrm{km}}{\mathrm{s}\cdot\mathrm{Mpc}}$$* Estimate mean and the uncertainty of the velocity of a galaxy which is measured to be $(500 \pm 100)\,\mathrm{Mpc}$ awayusing monte carlo methods<jupyter_code># %load 04_01_numpy_solutions/exercise_hubble.py<jupyter_output><empty_output><jupyter_text>Simple io functions<jupyter_code>idx = np.arange(100)
x = np.random.normal(0, 1e5, 100)
y = np.random.normal(0, 1, 100)
n = np.random.poisson(20, 100)
idx.shape, x.shape, y.shape, n.shape
np.savetxt(
    'data.txt',
    np.column_stack([idx, x, y, n]),
)
!head data.txt
# Load back the data, unpack=True is needed to read the data columnwise and not row-wise
idx, x, y, n = np.genfromtxt('data.txt', unpack=True)

idx.dtype, x.dtype<jupyter_output><empty_output><jupyter_text>Problems* Everything is a float* Way larger file than necessary because of too much digits for floats* No column names  Numpy recarrays* Numpy recarrays can store columns of different types* Rows are addressed by integer index* Columns are addressed by stringsSolution for our io problem → Column names, different types<jupyter_code># for more options on formatting see
# https://docs.scipy.org/doc/numpy/reference/generated/numpy.savetxt.html

data = np.savetxt(
    'data.csv',
    np.column_stack([idx, x, y, n]),
    delimiter=',', # true csv file
    header=','.join(['idx', 'x', 'y', 'n']),
    fmt=['%d', '%.4g', '%.4g', '%d'],  # One formatter for each column
)
!head data.csv
data = np.genfromtxt(
    'data.csv',
    names=True, # load column names from first row
    dtype=None, # Automagically determince best data type for each column
    delimiter=',',
)
data[:10]
data[0]
data['n']
data.dtype<jupyter_output><empty_output><jupyter_text>Linear algebraNumpy offers a lot of linear algebra functionality, mostly wrapping LAPACK<jupyter_code># symmetric matrix, use eigh
# If not symmetric, use eig
mat = np.array([
    [4, 2, 0],
    [2, 1, -3],
    [0, -3, 4]
])

eig_vals, eig_vecs = np.linalg.eigh(mat)

eig_vals, eig_vecs
np.linalg.inv(mat)<jupyter_output><empty_output><jupyter_text>Numpy matricesNumpy also has a matrix class, with operator overloading suited for matrices<jupyter_code>mat = np.matrix(mat)
mat.T
mat ** 2
mat * 5
mat.I
mat * np.matrix([1, 2, 3]).T<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.13877007888983936, 'fraction_numerical': 0.03526280771789754, 'mean_word_length': 4.995441595441595, 'pattern_counts': {'":': 0, '<': 129, '<?xml version=': 0, '>': 129, 'https://': 2, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 1, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '32', 'dir': 'jupyter-structured-clean-dedup', 'id': '465488', 'license': 'MIT', 'n_tokens_mistral': 4165, 'n_tokens_neox': 3891, 'n_words': 1085, 'path': 'notebooks/04_01_numpy.ipynb', 'repo_name': 'Python4AstronomersAndParticlePhysicists/PythonWorkshop-ICE', 'size': '10521'}