<jupyter_start><jupyter_text>Advanced NumPy usage<jupyter_code>from random import randint

gold = [randint(0, 100) for i in range(10000000)]
%%timeit
silver = []
for gold_value in gold:
    silver.append(12 * gold_value)
%timeit silver = [12 * gold_value for gold_value in gold]
import numpy as np
gold = np.random.randint(0, 100, (10000000,))
%%timeit
silver = np.empty(gold.shape)
for i in range(gold.shape[0]):
    silver[i] = 12 * gold[i]
%timeit silver = 12 * gold
loots = np.random.randint(0, 100, (10,))
shares = np.array([.7, 1., 1.3, 2., 1.])
loots * shares
gold = np.random.randint(0, 100, (10,))
gold
ratio = np.array(12)
ratio_wide = ratio * np.ones(gold.shape)
ratio_wide
gold * ratio_wide
loots
loots.shape
loots[:, np.newaxis]
loots[:, np.newaxis].shape
loots[:, None]
np.newaxis is None
loots = np.random.randint(0, 100, (10,))
shares = np.array([.7, 1., 1.3, 2., 1.])
loots[:, np.newaxis] * shares
names = np.array(['bard', 'mage', 'leader', 'paladin', 'ranger'])
# character attributes = health, speed, strength, mana
attributes = np.array([[72, 15, 48, 22],
                      [68, 8, 44, 87],
                      [65, 12, 56, 28],
                      [82, 7, 78, 0],
                      [67, 11, 62, 14]])
attributes
sorter = attributes[:, 0].argsort()
sorter
attributes[sorter]
names[sorter]
contributions = np.array([.7, 1., 1.3, 2., 1.])
attributes.T * contributions
loots = np.random.randint(-100, 100, (10,))
loots
loots > 0
loots[loots > 0]
loots[np.logical_and(loots >= 5, loots <= 30)]
dungeon_difficulties = np.array([[  7,   7,   5,   1,   5,   6,   9,   6,  10,   4,   2,   3],
                                 [ 10,  16,  12,  12,  18,   0,  16,  20,  12,   4,  12,  12],
                                 [  3,   3,  12,  15,   6,  12,  27,   6,   3,  30,  15,  27],
                                 [ 24,  12,  40,  32,  40,  12,  12,  16,  24,   4,  16,   4],
                                 [ 25,   0,  35,   0,  40,  50,  40,  20,  35,  25,  25,  35],
                                 [ 30,  24,  36,  48,  42,   6,  18,  54,  36,  54,  12,  24],
                                 [ 49,  28,  42,  70,  21,  42,   7,  42,  49,  42,   0,  63],
                                 [  8,  16,  72,  16,  24,  40,  40,  40,  64,  24,  24,   8],
                                 [ 81,  62,  18,  81,  81,  27,  27,  45,  18,  27,  90,  36],
                                 [ 80,  90,  20, 100,  80,  90,  50,  80,  70,  90,  60,  40]])
dungeon_difficulties[[1, 4, 7], [0, 2, 3]]
dungeon_difficulties[np.ix_([1, 4, 7], [0, 2, 3])]
np.ones(dungeon_difficulties.shape).astype(int)
import numpy as np
objects_in_dungeons = np.array([[1, 1, 1, 2, 1, 3, 1, 5, 7, 1, 1, 9],
                                [1, 1, 1, 1, 2, 2, 2, 1, 1, 1, 1, 10],
                                [1, 1, 1, 5, 1, 1, 3, 1, 2, 1, 1, 8],
                                [1, 1, 1, 1, 10, 1, 1, 2, 2, 2, 1, 13],
                                [1, 1, 1, 2, 3, 3, 5, 8, 1, 2, 4, 10],
                                [2, 3, 2, 2, 1, 5, 6, 1, 9, 2, 7, 6],
                                [2, 2, 2, 1, 8, 2, 2, 1, 6, 1, 5, 12],
                                [3, 2, 3, 5, 1, 1, 4, 2, 1, 7, 1, 9],
                                [3, 3, 2, 2, 4, 3, 6, 2, 8, 1, 1, 17],
                                [2, 2, 2, 3, 3, 1, 8, 1, 1, 1, 1, 18]])
np.isin(objects_in_dungeons, 8).nonzero()
objects_in_dungeons[np.isin(objects_in_dungeons, np.arange(10, 20))]
dungeon_list = np.array([1, 4, 7])
floor_list = np.array([0, 2, 3])
np.isin(dungeon_difficulties, dungeon_list)
%timeit dungeon_difficulties == 81
%timeit np.isin(dungeon_difficulties, np.array([81]))
dungeon_difficulties[dungeon_difficulties]
N = 10000
large_dungeons = (np.arange(10, 110, N)[:, None] * np.round(np.random.rand(N, 12), decimals=1)).astype(int)
import numpy as np
loots_per_character_per_floor = np.random.randint(0, 100, (10, 5, 12))
loots_per_character_per_floor[0]
loots_per_character_per_floor[:, :, -1]
loots_per_character_per_floor[..., -1]<jupyter_output><empty_output><jupyter_text>Indexing<jupyter_code>import numpy as np
np.set_printoptions(suppress=True)
devises = np.array([20, 1.8, .5, 12])
loots = np.concatenate([np.random.randint(0, 4, (10,))[:, None],
                        np.random.randint(0, 100, (10,))[:, None]],
                       axis=1)
loots
loots[:, 1][:, None] * devises
(loots[:, 1][:, None] * devises)[np.arange(len(loots)), loots[:, 0]]<jupyter_output><empty_output><jupyter_text>C-order vs F-order<jupyter_code>a = np.random.rand(100, 1000000)
%timeit  a.sum(axis=0)
b = a.astype(dtype=a.dtype, order='F')
%timeit  b.sum(axis=0)<jupyter_output>60.8 ms ± 70.7 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)<jupyter_text>Einsum<jupyter_code># character attributes = health, speed, strength, mana
attributes = np.array([[72, 15, 48, 22],
                      [68, 8, 44, 87],
                      [65, 12, 56, 28],
                      [82, 7, 78, 1],
                      [67, 11, 62, 14]])
attributes.shape
loots = np.random.randint(0, 100, (3, 5))
loots
loots[:, :, None] / attributes
np.einsum('ij, jk -> ijk', loots, 1 / attributes)<jupyter_output><empty_output><jupyter_text>View vs copy<jupyter_code>a = np.arange(18).reshape(3, 6)
a
a[np.arange(3), [1, 2, 3]] = -1
a
a[np.mod(a, 3) == 0] = -1
a
np.arange(len(loots)), loots[:, 0]
import numpy as np
a = np.arange(27).reshape((9, 3))
a<jupyter_output><empty_output><jupyter_text>Under the hood<jupyter_code>vectorization -> splitting loot part depending on character weight / contribution<jupyter_output><empty_output><jupyter_text>Vectorization The maze of indexing<jupyter_code>normalization in one line<jupyter_output><empty_output><jupyter_text>Broadcasting<jupyter_code>Crazy things with indexing
using weights / weighting samples<jupyter_output><empty_output><jupyter_text>Ufuncs<jupyter_code>import numpy as np
# closest value from a reference list
ref_list = np.array([2, 5, 12])
def closest(a, ref):
    return ref[(np.abs(ref - a)).argmin()]    
x = np.random.randint(0, 20, size=(1000))
%timeit np.apply_along_axis(closest, 1, x[:, None], ref_list)
closest_vect = np.vectorize(closest, excluded=['ref'])
%timeit closest_vect(x, ref=ref_list)
def closest_numpy(arr, ref):
    return ref[(np.abs(arr[:, None] - ref)).argmin(axis=1)]
ex = np.arange(5)
np.abs(ex[:, None] - ref_list)
%timeit closest_numpy(x, ref=ref_list)
from numba import vectorize, int64
@vectorize([int64(int64)], nopython=True)
def closest_numba(a):
    ref = ref_list
    return ref[(np.abs(ref - a)).argmin()]    
%timeit closest_numba(x)
ex = np.arange(15)
closest_numba.at(ex, ex < 10)
ex
np.add.reduce(ex)
ex.sum()<jupyter_output><empty_output><jupyter_text>Numba and Numexpr<jupyter_code>import numpy as np
np.set_printoptions(suppress=True)<jupyter_output><empty_output><jupyter_text>$ d_{euclid}(a, b) = \sqrt{||a-b||^2} =  \sqrt{ (a_1 - b_1)^2 + (a_2 - b_2)^2}  $<jupyter_code>N = 1000
A = np.random.rand(N, 2) + 1
B = np.random.rand(N, 2) + 1
def euclid(a, b):
    return np.sqrt(np.dot(a - b, (a - b).T))
def naive_euclid(X, Y):
    res = np.empty((len(X), len(Y)))
    for i in range(len(X)):
        for j in range(len(Y)):
            res[i, j] = euclid(X[i], Y[j])
    return res
%timeit naive_euclid(A, B)
def vectorized_euclid(X, Y):
    diff = X[np.newaxis] - Y[:, np.newaxis]
    return np.sqrt(np.einsum('ijk, ijk -> ji', diff, diff))
np.allclose(naive_euclid(A, B), vectorized_euclid(A, B))
%timeit vectorized_euclid(A, B)
from sklearn.metrics import pairwise_distances
np.allclose(pairwise_distances(A, B), vectorized_euclid(A, B))
%timeit pairwise_distances(A, B, metric='euclidean')
from numba import jit
@jit
def numba_euclid(X, Y):
    res = np.empty((len(X), len(Y)))
    for i in range(len(X)):
        for j in range(len(Y)):
            res[i, j] = np.sqrt(np.dot(X[i] - Y[j], X[i] - Y[j]))
    return res
%timeit numba_euclid(A, B)
@jit
def vectorized_euclid_numba(X, Y):
    diff = X[np.newaxis] - Y[:, np.newaxis]
    return np.sqrt(np.einsum('ijk, ijk -> ji', diff, diff))
%timeit vectorized_euclid_numba(A, B)<jupyter_output>22.2 ms ± 80 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)<jupyter_text>$d_{poinc}(a, b) = \mbox{arccosh} \left(1 + 2 \frac{||a - b ||^2}{(1 - ||a||^2) (1 - ||b|| ^2 )} \right)$<jupyter_code>def poinc(a, b):
    return np.arccosh(1+ 2* np.dot(a-b, a-b)/ ((1-np.dot(a,a))*(1-np.dot(b,b))))
def naive_poinc(X, Y):
    res = np.empty((len(X), len(Y)))
    for i in range(len(X)):
        for j in range(len(Y)):
            res[i, j] = poinc(X[i], Y[j])
    return res
%timeit naive_poinc(A, B)
def numpy_poinc(X, Y):
    diff = X[np.newaxis] - Y[:, np.newaxis]
    num = np.einsum('ijk, ijk -> ji', diff, diff)
    den1 = 1 - np.einsum('ij, ji->i', X, X.T)
    den2 = 1 - np.einsum('ij, ji->i', Y, Y.T)
    den = den1[:, np.newaxis] * den2[np.newaxis]
    return np.arccosh(1 + 2 * num / den)    
%timeit numpy_poinc(A, B)
@jit(nopython=True)
def numba_poinc(X, Y):
    def poinc(a, b):
        return np.arccosh(1+ 2* np.dot(a-b, a-b)/ ((1-np.dot(a,a))*(1-np.dot(b,b))))
    res = np.empty((len(X), len(Y)))
    i = 0
    for i in range(len(X)):
        for j in range(len(Y)):
            res[i, j] = poinc(X[i], Y[j])
    return res
%timeit numba_poinc(A, B)<jupyter_output>220 ms ± 7.53 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)<jupyter_text>Strides  Windowed average<jupyter_code>speed = np.random.rand(300000) * 75
def windowed_average(a, n):
    av = np.empty(len(a) - n + 1)
    for i in range(len(a) - n):
        av[i] = a[i:i+n].sum() / n
    return av
windowed_average(np.arange(30), 5)
%timeit windowed_average(speed, 5)
%timeit windowed_average_conv(speed, 5)
%timeit windowed_average_cumsum(speed, 5)
def windowed_average_conv(a, n):
    return np.convolve(a, np.ones((n,))/n, mode='valid')
windowed_average_conv(np.arange(30), 5)
def windowed_average_cumsum(a, n):
    cumsum = np.cumsum(np.insert(a, 0, 0)) 
    return (cumsum[n:] - cumsum[:-n]) / float(n)
windowed_average_cumsum(np.arange(30), 5)
from numpy.lib.stride_tricks import as_strided

def windowed_median_strides(a, n):
    stride = a.strides[0]
    return np.median(as_strided(a, shape=[len(a) - n + 1, n], strides=[stride, stride]), axis=1)
windowed_median_strides(np.arange(30), 5)
%timeit windowed_median_strides(speed, 5)
from numba import guvectorize, float64, int64
@guvectorize(['void(float64[:], int64, float64[:])'], '(n),()->(n)', nopython=True)
def windowed_median(a, n, out):
    for i in range(len(a) - n + 1):
        out[i] = np.median(a[i: i + n])
windowed_median(np.arange(30), 5)
%timeit windowed_median(speed, 5)
@guvectorize(['void(float64[:], int64, float64[:])'], '(n),()->(n)', nopython=True)
def windowed_average_numba(a, window_width, out):
    asum = 0.0
    count = 0
    for i in range(window_width):
        asum += a[i]
        count += 1
        out[i] = asum / count
    for i in range(window_width, len(a)):
        asum += a[i] - a[i - window_width]
        out[i] = asum / count
windowed_average_numba(np.arange(30), 5)
%timeit windowed_average_numba(speed, 5)<jupyter_output>421 µs ± 5.76 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
{'fraction_non_alphanumeric': 0.17915513083861748, 'fraction_numerical': 0.06850049120300081, 'mean_word_length': 2.926367461430575, 'pattern_counts': {'":': 0, '<': 51, '<?xml version=': 0, '>': 61, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 2, 'word_list_counts': {'cursed_substrings.json': 36, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '13', 'dir': 'jupyter-structured-clean-dedup', 'id': '513384', 'license': 'MIT', 'n_tokens_mistral': 5219, 'n_tokens_neox': 4676, 'n_words': 1138, 'path': 'NumPy/code_samples.ipynb', 'repo_name': 'ewjoachim/slides', 'size': '11197'}