The NumPy array¶
A structure for efficient numerical computation¶
Stéfan van der Walt (@stefanvdwalt)
University of California, Berkeley
Advanced Scientific Programming in Python
Summer School 2016 held in Reading, UK
Accompanying problem sets are on the Wiki at: https://python.g-node.org/wiki/numpy
Num-what?¶
This talk discusses some of the more advanced NumPy features. If you’ve never seen NumPy before, you may have more fun doing the tutorial at
http://mentat.za.net/numpy/intro/intro.html
or studying the course notes at
http://scipy-lectures.github.io
You can always catch up by reading:
The NumPy array: a structure for efficient numerical computation. Stéfan van der Walt, S. Chris Colbert and Gaël Varoquaux. In IEEE Computing in Science Engineering, March/April 2011.
Nicolas Rougier also has a beautiful tutorial at
Setup¶
import numpy as np # We always use this convention,
# also in the documentation
print(np.__version__) # version 1.9 or greater
ATLAS is a fast implementation of BLAS (Basic Linear Algebra Routines). On
OSX you have Accelerate; students can get Intel's MKL for free. On Ubuntu,
install libatlas3-base.
Make use of IPython's powerful features! TAB-completion, documentation, source inspection, timing, cpaste, etc.
print(np.show_config()) # got ATLAS / Accelerate / MKL?
Taking a look at numpy/core/include/numpy/ndarraytypes.h:
ndarray typedef struct PyArrayObject {
PyObject_HEAD
char *data; /* pointer to data buffer */
int nd; /* number of dimensions */
npy_intp *dimensions; /* size in each dimension */
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 struct */
int flags; /* Flags */
PyObject *weakreflist; /* For weakreferences */
} PyArrayObject ;
A homogeneous container¶
char *data; /* pointer to data buffer */
Data is just a pointer to bytes in memory:
x = np.array([1, 2, 3])
x.dtype
x.ctypes.data
import ctypes
list(ctypes.string_at(x.ctypes.data, x.size * x.itemsize))
Dimensions¶
int nd;
npy_intp *dimensions; /* number of dimensions */
/* size in each dimension */
x = np.array([])
np.array(0).shape
x = np.random.random((3, 2, 3, 3))
x.shape
x.ndim
Data type descriptors¶
PyArray_Descr * descr ; /* Pointer to type struct */
Common types in include int, float, bool:
np.array ([-1, 0, 1], dtype=int)
np.array([-1, 0, 1], dtype=float)
np.array([-1 , 0, 1], dtype=bool)
Each item in the array has to have the same type (occupy a fixed nr of bytes in memory), but that does not mean a type has to consist of a single item:
dt = np.dtype([('value', np.int), ('status', np.bool)])
np.array([(0, True), (1, False)], dtype=dt)
This is called a structured array.
Strides¶
npy_intp *strides; /* bytes to jump to get */
/* to the next element */
x = np.arange(12).reshape((3, 4))
x
x.dtype
x.dtype.itemsize
x.strides
We see (4 * itemsize, itemsize) or (skip_bytes_row, skip_bytes_col).
Example: Taking the transpose¶
We can take the transpose of the array without copying any data by simply manipulating shape and strides. What should they be?
Flags¶
int flags; /* Flags */
x = np.array([1, 2, 3])
z = x[::2]
x.flags
z.flags
Structured arrays¶
Like we saw earlier, each item in an array has the same type, but that does not mean a type has to consist of a single item:
dt = np.dtype([('value', np.int), ('status', np.bool)])
np.array([(0, True), (1, False)], dtype=dt)
This is called a structured array, and is accessed like a dictionary:
x = np.array([(0, True), (1, False)], dtype=dt)
x['value']
x['status']
Reading structured data from file¶
Reading this kind of data can be troublesome:
while ((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:
dt = 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:
data['position']['az']
Broadcasting¶
Combining of differently shaped arrays without creating large intermediate arrays:
x = np.arange(4)
print(x)
print(x + 3)
See docstring of np.doc.broadcasting for more detail.
Broadcasting (2D)¶
a = np.arange(12).reshape((3, 4))
b = np.array([1, 2, 3])[:, np.newaxis]
print(a)
print()
print(b)
a + b
Broadcasting (3D)¶
x = np.zeros((3, 5))
y = np.zeros(8)
z = x[..., np.newaxis]
x = np.zeros((3, 5))
y = np.zeros(8)
z = x[..., np.newaxis]
X shape = (3, 5, 1)
Y shape = ( 8)
Z shape = (3, 5, 8)
Broadcasting rules¶
The broadcasting rules are straightforward: compare dimensions, starting from
the last. Match when either dimension is one or None, or if dimensions are
equal:
Scalar 2D 3D Bad
( ,) (3, 4) (3, 5, 1) (3, 5, 2)
(3,) (3, 1) ( 8) ( 8)
---- ------ --------- ---------
(3,) (3, 4) (3, 5, 8) XXXExplicit broadcasting¶
x = np.zeros((3, 5, 1))
y = np.zeros((3, 5, 8))
xx, yy = np.broadcast_arrays(x, y)
xx.shape
np.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]]np.broadcast_to(x, (3,5,8)).shape
Fancy indexing¶
An ndarray can be indexed in two ways:
- Using slices and scalars
- Using other ndarrays ("fancy indexing")
x = np.arange(9).reshape((3, 3))
x
x[:, [1, 1, 2]]
Which is equivalent to:
np.array((x[:, 1], x[:, 1], x[:, 2])).T
Output shape of an indexing op¶
- Broadcast all index arrays against one another.
- Use the dimensions of slices as-is.
Start with the following input:
x = np.array([[0, 1, 2], [3, 4, 5], [6, 7, 8]])
x.shape
What would happen if we do
idx0 = np.array([[0, 1],
[1, 2]]) # row indices
idx1 = np.array([[0, 1]]) # column indices
x[idx0, idx1]
Output shape of indexing op (cont'd)¶
idx0.shape, idx1.shape
a, b = np.broadcast_arrays(idx0, idx1)
a
b
x
x[idx0, idx1]
Output shape of an indexing op (cont'd)¶
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)
Predict the output shape of
x[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)¶
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)
Broadcast index1 against index2:
(3, 2) # shape of index1
(1, 2) # shape of index2
------
(3, 2)The shape of x[5:10, index_one, :, index_two] is
(3, 2, 5, 16)
Jack's dilemma¶
Indexing and broadcasting are intertwined, as we’ll see in the following example. One of my favourites from the NumPy mailing list:
email
Date: Wed, 16 Jul 2008 16:45:37 -0500
From: Jack Cook
To: <numpy-discussion@scipy.org>
Subject: Numpy Advanced Indexing QuestionGreetings,
I have an I,J,K 3D volume of amplitude values at regularly sampled time 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 a subvolume at a constant +/- K window around the slice. Is there an easy way to do this using advanced indexing or some other method? Thanks in advanced for your help.
-- Jack
Test setup¶
ni, nj, nk = (10, 15, 20) # dimensions
half_width = 3 # take 7 elements
Make a fake data block such that block[i, j, k] == k for all i, j, k.
block = np.empty((ni, nj, nk) , dtype=int)
block[:] = np.arange(nk)[np.newaxis, np.newaxis, :]
Pick out a random fake horizon in k:
rng = np.random.RandomState(0)
k = rng.randint(5, 15, size=(ni, nj))
k
Solution¶
These two indices assure that we take a slice at each (i, j) position
idx_i = np.arange(ni)[:, np.newaxis, np.newaxis]
idx_j = np.arange(nj)[np.newaxis, :, np.newaxis]
This is the substantitive part that picks out the window:
idx_k = k[:, :, np.newaxis] + np.arange(-half_width, half_width + 1)
slices = block[idx_i, idx_j, idx_k] # Slice
Apply the broadcasting rules:
(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¶
slices = block[idx_i, idx_j, idx_k]
slices.shape
(10, 15, 7)
Now verify that our window is centered on k everywhere:
slices[:, :, 3]
np.all(slices[:, :, 3] == k)
The __array_interface__¶
Any object that exposes a suitable dictionary named
__array_interface__ may be converted to a NumPy array. This is
handy for exchanging data with external libraries. The array interface
has the following important keys (see
http://docs.scipy.org/doc/numpy/reference/arrays.interface):
- shape
- typestr
- data: (20495857, True); 2-tuple—pointer to data and boolean to indicate whether memory is read-only
- strides
- version: 3
Stride Tricks: Repeating sequence¶
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
Did you know?¶
NumPy implements a Diophantine (integer) equation solver for np.shares_memory to calculate whether two arrays possibly overlap in memory.
sum(a[i] * x[i] for i in range(n)) == b
a[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?
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)
Questions, discussion and exercises¶
%reload_ext load_style
%load_style css/notebook.css