<jupyter_start><jupyter_text>Numpy laboratory<jupyter_code>import numpy as np

import matplotlib.pyplot as plt
%matplotlib inline<jupyter_output><empty_output><jupyter_text>1. Numpy basic exercisesUse vectorization and avoid for loops in all exercises. Implement each exercise as a single function and write test cases for the function. 1.1 Implement standardization for 2D arrays.\begin{equation*}X_{std} = \frac{X - \mu}{\sigma},\end{equation*}where $\mu$ is the mean and $\rho$ is the standard deviation of the array elements.<jupyter_code>def standardize(X):
    return (X - X.mean(axis=0)) / X.std(axis=0)

x_std = standardize(np.array([[1, 2], [2, 3], [1, 2]]))
x_std<jupyter_output><empty_output><jupyter_text>1.2 Implement normalization for 2D arrays.\begin{equation*}X_{norm} = \frac{X - X_{min}}{X_{max} - X_{min}}\end{equation*}<jupyter_code>def normalize(X):
    xmin = X.min(axis=0)
    xmax = X.max(axis=0)
    return (X - xmin) / (xmax - xmin)

X = np.arange(6).reshape(2, 3)
X[1, 2] = -5
print(X)
normalize(X)<jupyter_output>[[ 0  1  2]
 [ 3  4 -5]]<jupyter_text>1.3 Implement the softmax function.$$x_i \mapsto \frac{\exp(x_i)}{\sum_{j=1}^n \exp(x_j)}$$<jupyter_code>def softmax(X):
    # compute softmax along the last dimension
    expX = np.exp(X)
    return expX / expX.sum(axis=-1)[:, None]

X = np.arange(12).reshape(3, -1)

softmax(X)<jupyter_output><empty_output><jupyter_text>2. VectorizationRewrite the following examples into vectorized solutions (no for loops and list comprehensions). 2.1 Row-wise Euclidean normWrite a function which has one parameter, a 2D array and it returns a vector of row-wise Euclidean norms of the input. Use `numpy` operations and vectorization, avoid `for` loops. The solution below is a _bad_ solution.<jupyter_code>def rowwise_norm(X):
    def my_dot(x, y):
        result = 0.0
        for i in range(len(x)):
            result += x[i] * y[i]
        return result
    return np.array([np.sqrt(my_dot(x, x)) for x in X])

X = np.arange(5)[:, None]*np.ones((5, 3));
print(X)
print(rowwise_norm(X))
print(rowwise_norm([[1], [-1], [1], [-1]]))<jupyter_output>[[0. 0. 0.]
 [1. 1. 1.]
 [2. 2. 2.]
 [3. 3. 3.]
 [4. 4. 4.]]
[0.         1.73205081 3.46410162 5.19615242 6.92820323]
[1. 1. 1. 1.]<jupyter_text>Solution<jupyter_code>def rowwise_norm(X):
    X = np.array(X)
    return np.sqrt((X * X).sum(axis=1))

X = np.arange(5)[:, None] * np.ones((5, 3));
print(X)
print(rowwise_norm(X))
print(rowwise_norm([[1], [-1], [1], [-1]]))<jupyter_output>[[0. 0. 0.]
 [1. 1. 1.]
 [2. 2. 2.]
 [3. 3. 3.]
 [4. 4. 4.]]
[0.         1.73205081 3.46410162 5.19615242 6.92820323]
[1. 1. 1. 1.]<jupyter_text>2.2 ChessboardWrite a function which has one parameter, a positive integer $n$, and returns an $n\times n$ array of $\pm1$ values like a chessboard: $M_{i,j} = (-1)^{i+j}$.<jupyter_code>def chessboard(n):
    return np.array([[(-1)**(i + j) for j in range(n)] for i in range(n)])

chessboard(5)<jupyter_output><empty_output><jupyter_text>Solution<jupyter_code>def chessboard(n):
    return (-1) ** np.arange(n * n).reshape(n, n)

chessboard(3)<jupyter_output><empty_output><jupyter_text>3. Broadcast quizDo the following operations work and if so, what is the shape of the resulting array?* Try to figure it out before evaluating the cell.<jupyter_code>np.ones(3) + np.ones((3,3))
np.ones(3) + np.ones((4, 3))
# np.ones(3) + np.ones((3, 4))  # fails
np.ones(3)[:, None] + np.ones((3, 4))
np.ones((1, 2, 3)) + np.ones((1, 3))[:, None, :]
np.ones((1, 2, 3)) + np.ones((1, 3))<jupyter_output><empty_output><jupyter_text>Can you broadcast these arrays such that they can be added? What will be the shape of the result?<jupyter_code># np.ones((4,20,3)) + np.ones((5, 3))  # fails<jupyter_output><empty_output><jupyter_text>4. Numpy advanced exercises 4.1 BlockmatrixWrite a function named __`blockmatrix`__ that produces the following (block) matrix:$$\left(\begin{array}{ccc|ccc} 1 & & 0& 0 & \cdots & 0 \\ & \ddots & & \vdots & & \vdots \\  0& & 1 & 0 & \cdots & 0 \\\hline  0 & \cdots & 0 & 1 & \cdots & 1 \\  \vdots & & \vdots & \vdots & & \vdots \\  0 & \cdots & 0 & 1 & \cdots & 1\end{array}\right)$$The function should have 2 positive integer parameters, the size of the first square block and the size of the last square block. The other two blocks should have the appropriate size (may be rectangle).The first block is an indentity matrix, the last is a constant $1$ matrix.Return the resulted matrix.Use matrix initializers: `ones`, `zeros`, `eye` and concatenation.<jupyter_code>def blockmatrix(n, m):
    return np.concatenate(
        (
            np.concatenate((np.eye(n), np.zeros((n, m))), axis=1),  # upper block
            np.concatenate((np.zeros((m, n)), np.eye(m)), axis=1),  # lower block
        ),
        axis=0
    )

blockmatrix(2, 5)<jupyter_output><empty_output><jupyter_text>4.2 Blockmatrix from arbitrary square matricesWrite a blockmatrix function that takes any number of square matrices and returns a blockmatrix with these matrices in the diagonal.<jupyter_code>def any_blockmatrix(*matrices):
    sum_width = sum(m.shape[1] for m in matrices)
    before = 0
    X = np.concatenate((matrices[0], np.zeros((matrices[0].shape[0], sum_width-matrices[0].shape[1]))), axis=1)
    before += matrices[0].shape[1]
    for m in matrices[1:]:
        after = sum_width - before - m.shape[1]
        height = m.shape[0]
        X = np.concatenate((
            X,
            np.concatenate((np.zeros((height, before)), m, np.zeros((height, after))), axis=1)
        ))
        before += m.shape[1]
    return X

A = np.arange(9).reshape(3, 3)
B = np.diag((2, 3))
any_blockmatrix(A, B, A, B)<jupyter_output><empty_output><jupyter_text>4.3 DerivativeWrite a function which numerically derivates a $\mathbb{R}\mapsto\mathbb{R}$ function. Use the forward finite difference.The input is a 1D array of function values, and optionally a 1D vector of abscissa values. If not provided then the abscissa values are unit steps.The result is a 1D array with the length of one less than the input array.Use `numpy` operations instead of `for` loop in contrast to the solution below.<jupyter_code>def derivate(f, x=None):
    if x is None:
        x = np.arange(len(f))
    return np.array([(f[i+1] - f[i]) / (x[i+1] - x[i]) for i in range(len(x) - 1)])

derivate(np.arange(10)**2)
x = np.arange(10)
plt.plot(x, x**2)
plt.plot(x[:-1], derivate(x**2, x))<jupyter_output><empty_output><jupyter_text>Solution<jupyter_code>def derivate(f, x=None):
    if x is None:
        x = np.arange(len(f))
    return np.diff(f) / np.diff(x)
x = np.arange(10)
plt.plot(x, x**2)
plt.plot(x[:-1], derivate(x**2, x))<jupyter_output><empty_output><jupyter_text>4.4 Birthday problemIn probability theory, the birthday problem or birthday paradox concerns the probability that, in a set of n randomly chosen people, some pair of them will have the same birthday. By the pigeonhole principle, the probability reaches 100% when the number of people reaches 367 (since there are only 366 possible birthdays, including February 29). However, 99.9% probability is reached with just 70 people, and 50% probability with 23 people. These conclusions are based on the assumption that each day of the year (excluding February 29) is equally probable for a birthday. -- [Wikipedia](https://en.wikipedia.org/wiki/Birthday_problem)Write a function that simulates this problem for variable $n$. Your function should take $n$ and an experiment count as its parameter and sample experiment count times and return the ratio of "birthday collisions" (how many times there were at least two birthdays on the same day).Run it for different $n$ values with at least 1000 experiment count and plot the results. You can add a grid and check for the 50% probability with:    fig, ax = plt.subplots()    ax.plot(x)    ax.grid()<jupyter_code>def simulate_birthday_problem(n, iter_no=1000):
    """Simulate the birthday problem for a group of N people
    Generates N birthdays with uniform probability and returns
    the proportion of iterations with birthday collisions.
    """
    collision = 0
    for _ in range(iter_no):
        birthdays = np.random.randint(1, 366, size=n)
        if np.unique(birthdays).shape[0] < n:  # collision
            collision += 1
    return collision / iter_no

max_n = 50
fig, ax = plt.subplots(figsize=(16, 8))
ax.set_yticks(np.arange(0, 1.0, .1))
ax.grid(True)

# number of simulations to run
iter_nums = [20, 100, 1000]

for iter_no in iter_nums:
    x = []
    for n in range(2, max_n):
        x.append(simulate_birthday_problem(n, iter_no=iter_no))
    ax.plot(x)
    
ax.legend(iter_nums)

# notice that the plot becomes less noisy as we run more simulations<jupyter_output><empty_output><jupyter_text>4.5 Horner's methodImplement the [Horner's method](https://en.wikipedia.org/wiki/Horner%27s_methodDescription_of_the_algorithm) for evaluating polynomials. The first input is a 1D array of numbers, the coefficients, from the constant coefficient to the highest order coefficent. The second input is the variable $x$ to subsitute. The function should work for all type of variables: numbers, arrays; the output should be the same type array as the input, containing the elementwise polynomial values.<jupyter_code>def horner(C, x):
    y = np.zeros_like(x)
    for c in C[::-1]:
        y = y*x + c
    return y

C = [2, 0, 1] # 2 + 0*x + 1*x^2
print(horner(C, 3)) # 2 + 3^2
print(horner(C, [3, 3]))<jupyter_output>11
[11 11]
StarcoderdataJupyterStructuredCleanDedup