<jupyter_start><jupyter_text>The NumPy array A structure for efficient numerical computation<NAME> (@stefanvdwalt)University of California, BerkeleyAdvanced Scientific Programming in PythonSummer School 2018 held in Melbourne, AustraliaAccompanying **problem sets** are linked from the Wiki at https://www.melbournebioinformatics.org.au/aspp-asia-pacific/numpy   Num-what?This talk discusses some of the more advanced NumPy features. Most of youhave some experience with NumPy, but if not you can always catch up byfollowing the tutorial athttp://mentat.za.net/numpy/intro/intro.htmlor studying the course notes athttp://scipy-lectures.orgIf you want to dig a little deeper, the following paper is available from Arxiv:>  **The NumPy array: a structure for efficient numerical computation**. Stéfan>  <NAME>, <NAME> and <NAME>. In IEEE Computing in>  Science Engineering, March/April 2011.<NAME> has a beautiful tutorial athttps://github.com/rougier/numpy-tutorial  Setup<jupyter_code>import numpy as np  # We always use this convention,
                    # also in the documentation
print(np.__version__)   # version 1.11 or greater<jupyter_output>1.13.3<jupyter_text>NumPy should be linked to a fast implementation of BLAS (Basic Linear Algebra Routines)such as OpenBLAS or ATLAS.  On Ubuntu, e.g., install ``libatlas3-base``.  If you installed NumPy via `pip` or `conda`, this should already be the case.<jupyter_code>print(np.show_config()) # got ATLAS / Accelerate / MKL?<jupyter_output>blas_mkl_info:
    libraries = ['mkl_rt', 'pthread']
    library_dirs = ['/home/student/anaconda3/lib']
    define_macros = [('SCIPY_MKL_H', None), ('HAVE_CBLAS', None)]
    include_dirs = ['/home/student/anaconda3/include']
blas_opt_info:
    libraries = ['mkl_rt', 'pthread']
    library_dirs = ['/home/student/anaconda3/lib']
    define_macros = [('SCIPY_MKL_H', None), ('HAVE_CBLAS', None)]
    include_dirs = ['/home/student/anaconda3/include']
lapack_mkl_info:
    libraries = ['mkl_rt', 'pthread']
    library_dirs = ['/home/student/anaconda3/lib']
    define_macros = [('SCIPY_MKL_H', None), ('HAVE_CBLAS', None)]
    include_dirs = ['/home/student/anaconda3/include']
lapack_opt_info:
    libraries = ['mkl_rt', 'pthread']
    library_dirs = ['/home/student/anaconda3/lib']
    define_macros = [('SCIPY_MKL_H', None), ('HAVE_CBLAS', None)]
    include_dirs = ['/home/student/anaconda3/include']
None<jupyter_text>Make use of Jupyter's powerful features!- TAB-completion- documentation (use `np.some_function?` or `some_func??` for source)- magics (e.g., `%timeit`, `%prun`, `%magic`)- cells magics (e.g., `%%bash`, `%%file`)  The NumPy ndarray Revision: Structure of an array   Taking a look at ``numpy/core/include/numpy/ndarraytypes.h``:```cppndarray typedef struct PyArrayObject {    PyObject_HEAD    char *data;             /* pointer to data buffer */    int nd;                 /* number of dimensions, also called `ndim` */    npy_intp *dimensions;   /* size in each dimension, also called `shape`  */    npy_intp *strides;      /* bytes to jump to get                             * to the next element in each dimension                             */    PyObject *base;         /* pointer to original array                            /* decref this object upon deletion. */    PyArray_Descr *descr;   /* pointer to type structure, called `dtype` */    int flags;              /* flags */    PyObject *weakreflist;  /* for weak references */} PyArrayObject ;```  A homogeneous container ```cppchar *data;  /* pointer to data buffer */``` Data is just a pointer to bytes in memory:<jupyter_code>x = np.array([1, 2, 3])
x.dtype
x.ctypes.data
import ctypes
list(ctypes.string_at(x.ctypes.data, x.size * x.itemsize))<jupyter_output><empty_output><jupyter_text>Dimensions ```cppint nd;npy_intp *dimensions; /* number of dimensions */                      /* size in each dimension */```<jupyter_code>x = np.array([])
np.array(0).shape
x = np.random.random((3, 2, 3, 3))
x.shape
x.ndim<jupyter_output><empty_output><jupyter_text>Data type descriptors ```cppPyArray_Descr * descr ;  /* Pointer to type struct */``` Common types in include ``int``, ``float``, ``bool``:<jupyter_code>np.array ([-1, 0, 1], dtype=int)
np.array([-1, 0, 1], dtype=float)
np.array([-1 , 0, 1], dtype=bool)<jupyter_output><empty_output><jupyter_text>Each item in the array has to have the same type (occupy a fixed nr of bytesin memory), but that does not mean a type has to consist of a single item:<jupyter_code>dt = np.dtype([('value', np.int), ('status', np.bool)])
np.array([(0, True), (1, False)], dtype=dt)<jupyter_output><empty_output><jupyter_text>This is called a structured array.  Strides ```cppnpy_intp *strides;  /* bytes to jump to get */                    /* to the next element */```<jupyter_code>x = np.arange(12).reshape((3, 4))
x
x.dtype
x.dtype.itemsize
x.strides<jupyter_output><empty_output><jupyter_text>We see ``(4 * itemsize, itemsize)`` or ``(skip_bytes_row, skip_bytes_col)``.   Example: Taking the transposeWe can take the transpose of the array *without* copying any data by simply manipulating ``shape`` and ``strides``.  What should they be?  Flags ```cppint flags;  /* Flags */```<jupyter_code>x = np.array([1, 2, 3])
z = x[::2]

x.flags
z.flags<jupyter_output><empty_output><jupyter_text>Structured arraysLike we saw earlier, each item in an array has the same type, but that does notmean a type has to consist of a single item:<jupyter_code>dt = np.dtype([('value', np.int), ('status', np.bool)])
np.array([(0, True), (1, False)], dtype=dt)<jupyter_output><empty_output><jupyter_text>This is called a structured array, and is accessed like a dictionary:<jupyter_code>x = np.array([(0, True), (1, False)], dtype=dt)
x['value']
x['status']<jupyter_output><empty_output><jupyter_text>Reading structured data from file Reading this kind of data can be troublesome: ```matlabwhile ((count > 0) && (n <= NumPoints))  % get time - I8 [ ms ]  [lw, count] = fread(fid, 1, 'uint32');  if (count > 0) % then carry on  uw = fread(fid, 1, 'int32');  t(1, n) = (lw + uw * 2^32) / 1000;  % get number of bytes of data  numbytes = fread(fid, 1, 'uint32');  % read sMEASUREMENTPOSITIONINFO (11 bytes)  m(1, n) = fread(fid, 1, 'float32'); % az [ rad ]  m(2, n) = fread(fid, 1, 'float32'); % el [ rad ]  m(3, n) = fread(fid, 1, 'uint8');   % region type  m(4, n) = fread(fid, 1, 'uint16');  % region ID  g(1, n) = fread(fid, 1, 'uint8');  numsamples = ( numbytes -12)/2; % 2 byte int  ...``` Using structured arrays:```pythondt = np.dtype([('time', np.uint64),               ('size', np.uint32),               ('position', [('az', np.float32),                             ('el', np.float32),                             ('region_type', np.uint8),                             ('region_ID', np.uint16)]),               ('gain', np.uint8),               ('samples', (np.int16, 2048))])data = np.fromfile(f, dtype=dt)``` We can then access this structured array as before:```pythondata['position']['az']``` If this is the type of data you work with most frequently,you should also look at Pandas (see [my DS tutorial](https://github.com/stefanv/ds_intro)).  BroadcastingCombining of differently shaped arrays without creating large intermediatearrays:<jupyter_code>x = np.arange(4)
print(x)
print(x + 3)<jupyter_output>[0 1 2 3]
[3 4 5 6]<jupyter_text>![Broadcasting in 1D](images/broadcast_1D.png)See docstring of ``np.doc.broadcasting`` for more detail.  Broadcasting (2D)<jupyter_code>a = np.arange(12).reshape((3, 4))
b = np.array([1, 2, 3])[:, np.newaxis]

print(a)
print()
print(b)
a + b<jupyter_output><empty_output><jupyter_text>Broadcasting (3D)<jupyter_code>x = np.zeros((3, 5))
y = np.zeros(8)
z = x[..., None] + y<jupyter_output><empty_output><jupyter_text>```X shape = (3, 5, 1)Y shape = (      8)-------------------Z shape = (3, 5, 8)```   Broadcasting rulesThe broadcasting rule is: compare dimensions, starting fromthe last. Match when either dimension is one or ``None``, or if dimensions areequal. ```Scalar    2D           3D           Bad( ,)     (3, 4)     (3, 5, 1)    (3, 5, 2)(3,)     (3, 1)     (      8)    (      8)----     ------     ---------    ---------(3,)     (3, 4)     (3, 5, 8)       XXX```  Explicit broadcasting<jupyter_code>x = np.zeros((3, 5, 1))
y = np.zeros((3, 5, 8))

xx, yy = np.broadcast_arrays(x, y)

xx.shape<jupyter_output><empty_output><jupyter_text>This should familiar to those of you who know "meshgrid" from MATLAB. ```pythonnp.broadcast_arrays([1, 2, 3],                    [[1], [2], [3]])``````[[1 2 3]   [[1 1 1] [1 2 3]    [2 2 2] [1 2 3]]   [3 3 3]]```<jupyter_code># Alternative API for broadcasting explicitly:

np.broadcast_to(x, (3,5,8));<jupyter_output><empty_output><jupyter_text>Note: when broadcasting, NumPy does *not* generate big new matrices.  See section later on stride tricks, and consider the effect of `stride = 0`.  Fancy indexingAn ndarray can be indexed in two ways: - Using slices and scalars - Using other ndarrays ("fancy indexing")<jupyter_code>x = np.arange(9).reshape((3, 3))
x
x[:, 2]  # standard, slice-scalar
x[2, 2]  # standard, scalar-scalar
x[:, [1, 1, 2]]  # fancy: slice-array<jupyter_output><empty_output><jupyter_text>Output shape of an indexing op Consider the following input:<jupyter_code>x = np.array([[0, 1, 2], [3, 4, 5], [6, 7, 8]])
x.shape<jupyter_output><empty_output><jupyter_text>What would happen if we do<jupyter_code>idx0 = np.array([[0, 1],
                 [1, 2]])  # row indices

idx1 = np.array([[0, 1]])  # column indices<jupyter_output><empty_output><jupyter_text>```pythonx[idx0, idx1]```  Fancy indexing rule1. Broadcast all index arrays against one another.2. Use the dimensions of slices as-is.  Output shape of indexing op (cont'd)<jupyter_code># x[idx0, idx1]
idx0.shape, idx1.shape
a, b = np.broadcast_arrays(idx0, idx1)
a
b
x
x[idx0, idx1]<jupyter_output><empty_output><jupyter_text>Output shape of an indexing op (cont'd)<jupyter_code>x = np.random.random((15, 12, 16, 3))

index_one = np.array([[0, 1], [2, 3], [4, 5]])
index_two = np.array([[0, 1]])
print(index_one.shape)
print(index_two.shape)<jupyter_output>(3, 2)
(1, 2)<jupyter_text>Predict the output shape of ```pythonx[5:10, index_one, :, index_two]``` Warning! When mixing slicing and fancy indexing, the*order* of the output dimensions is less easy to predict.Play it safe and **don't mix the two!**  Output shape of an indexing op (cont'd)<jupyter_code>x = np.random.random((15, 12, 16, 3))
index_one = np.array([[0, 1], [2, 3], [4, 5]])
index_two = np.array([[0, 1]])
print(index_one.shape)
print(index_two.shape)<jupyter_output>(3, 2)
(1, 2)<jupyter_text>Broadcast ``index1`` against ``index2``: ```(3, 2)   shape of index1(1, 2)   shape of index2------(3, 2)```<jupyter_code>x[5:10, index_one, :, index_two].shape<jupyter_output><empty_output><jupyter_text>Jack's dilemmaNow that we are able to predict the output of broadcasting and indexing operations, let's take a look at one of my favourite questions from the NumPy mailing list: ```emailDate: Wed, 16 Jul 2008 16:45:37 -0500From: <NAME>To: Subject: Numpy Advanced Indexing Question```Greetings,I have an I,J,K 3D volume of amplitude values at regularly sampledtime intervals. I have an I,J 2D slice which contains a time (K)value at each I, J location. What I would like to do is extract asubvolume at a constant +/- K window around the slice. Is there aneasy way to do this using advanced indexing or some other method?Thanks in advanced for your help.-- Jack   Test setup<jupyter_code>import numpy as np
ni, nj, nk = (10, 15, 20)  # dimensions
half_width = 3  # take 7 elements<jupyter_output><empty_output><jupyter_text>Make a fake data block such that ``block[i, j, k] == k`` for all i, j, k.<jupyter_code>block = np.empty((ni, nj, nk) , dtype=int)
block[:] = np.arange(nk)[np.newaxis, np.newaxis, :]
%matplotlib inline
import matplotlib.pyplot as plt

f, (ax0, ax1, ax2) = plt.subplots(1, 3, figsize=(20, 10))

ax0.imshow(block[:, :, 15]);
ax0.set_title('block[:, :, 15]')

ax1.imshow(block[:, 10, :])
ax1.set_title('block[:, 10, :]')

ax2.imshow(block[5, :, :])
ax2.set_title('block[5, :, :]');

plt.savefig('_slices.png', bbox_inches='tight')
plt.close()<jupyter_output><empty_output><jupyter_text>Pick out a random fake horizon in ``k``:<jupyter_code>rng = np.random.RandomState(0)

k = rng.randint(5, 15, size=(ni, nj))
k
print(block.shape)
print(half_width)<jupyter_output>(10, 15, 20)
3<jupyter_text>Pause to think: what shape should the result be?  SolutionThese two indices assure that we take a slice at each (i, j) position<jupyter_code>idx_i = np.arange(ni)[:, np.newaxis, np.newaxis]  # shape (ni, 1, 1)
idx_j = np.arange(nj)[np.newaxis, :, np.newaxis]  # shape (1, nj, 1)<jupyter_output><empty_output><jupyter_text>This is the substantitive part that picks out the window:<jupyter_code>idx_k = k[:, :, np.newaxis] + np.arange(-half_width, half_width + 1)
# What is the shape of idx_k?
slices = block[idx_i, idx_j, idx_k]  # Slice<jupyter_output><empty_output><jupyter_text>Apply the broadcasting rules:```pythonidx_i = np.arange(ni)[:, np.newaxis, np.newaxis]   shape (ni, 1, 1)idx_j = np.arange(nj)[np.newaxis, :, np.newaxis]   shape (1, nj, 1)idx_k = k[:, :, np.newaxis] + np.arange(-half_width, half_width + 1)slices[block_idx, idx_j, idx_k]``` ```(ni,  1, 1                 )   idx_i(1,  nj, 1                 )   idx_j(ni, nj, 2 * half_width + 1)   idx_k----------------------------(ni, nj, 7)  <-- this is what we wanted!```  Solution verification<jupyter_code>slices = block[idx_i, idx_j, idx_k]
slices.shape
(10, 15, 7)<jupyter_output><empty_output><jupyter_text>Now verify that our window is centered on ``k`` everywhere:<jupyter_code>slices[:, :, 3]
np.all(slices[:, :, 3] == k)<jupyter_output><empty_output><jupyter_text>The \_\_array\_interface\_\_Any object that exposes a suitable dictionary named``__array_interface__`` may be converted to a NumPy array. This ishandy for exchanging data with external libraries. The array interfacehas the following important keys (seehttp://docs.scipy.org/doc/numpy/reference/arrays.interface): - **shape** - **typestr** - **data**: (20495857, True); 2-tuple—pointer to data and boolean toindicate whether memory is read-only - **strides** - **version**: 3  Stride Tricks: Repeating sequence<jupyter_code>data = np.array([1, 2, 3, 4])
print(data.shape, data.strides)
N = 400000
repeat = np.lib.stride_tricks.as_strided(data, shape=(N, 4), strides=(0, 8))
repeat<jupyter_output><empty_output><jupyter_text>Side TriviaNumPy implements a Diophantine (integer) [equation solver](https://github.com/numpy/numpy/blob/master/numpy/core/src/private/mem_overlap.c) for `np.shares_memory` to calculate whether two arrays possibly overlap in memory.Q: Can you imagine why this (`np.shares_memory`) is a tricky problem to solve? ```sum(a * x[i] for i in range(n)) == ba[i] > 0, 0 <= x[i] <= ub[i]``` You can even access this solver from Python.  For example, given that a year has 365 days, how do we fit in months of 28, 30, and 31 days?<jupyter_code>from numpy.core.multiarray_tests import solve_diophantine

month_days = (28, 30, 31)
limits = (20, 20, 20)

solve_diophantine(month_days, limits, 365)
(2 * 28) + (1 * 30) + (9 * 31) 
solve_diophantine(month_days, (1, 20, 20), 365)
(1 * 28) + (4 * 30) + (7 * 31)<jupyter_output><empty_output><jupyter_text>Questions, discussion and exercises ---<jupyter_code>%reload_ext load_style
%load_style css/notebook.css<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.15522696455091367, 'fraction_numerical': 0.029508620133014786, 'mean_word_length': 4.277001703577513, 'pattern_counts': {'":': 0, '<': 133, '<?xml version=': 0, '>': 135, 'https://': 4, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '32', 'dir': 'jupyter-structured-clean-dedup', 'id': '592540', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 5864, 'n_tokens_neox': 5491, 'n_words': 1585, 'path': 'numpy_advanced.ipynb', 'repo_name': 'abbatux/aspp_melbourne_numpy', 'size': '15536'}