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
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
import numpy as np # We always use this convention,
# also in the documentation
print(np.__version__) # version 1.9 or greater
1.11.1
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?
openblas_lapack_info: library_dirs = ['/usr/local/lib'] language = c define_macros = [('HAVE_CBLAS', None)] libraries = ['openblas', 'openblas'] blas_opt_info: library_dirs = ['/usr/local/lib'] language = c define_macros = [('HAVE_CBLAS', None)] libraries = ['openblas', 'openblas'] blas_mkl_info: NOT AVAILABLE lapack_opt_info: library_dirs = ['/usr/local/lib'] language = c define_macros = [('HAVE_CBLAS', None)] libraries = ['openblas', 'openblas'] openblas_info: library_dirs = ['/usr/local/lib'] language = c define_macros = [('HAVE_CBLAS', None)] libraries = ['openblas', 'openblas'] None
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 ;
char *data; /* pointer to data buffer */
Data is just a pointer to bytes in memory:
x = np.array([1, 2, 3])
x.dtype
dtype('int64')
x.ctypes.data
18269632
import ctypes
list(ctypes.string_at(x.ctypes.data, x.size * x.itemsize))
[1, 0, 0, 0, 0, 0, 0, 0, 2, 0, 0, 0, 0, 0, 0, 0, 3, 0, 0, 0, 0, 0, 0, 0]
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
(3, 2, 3, 3)
x.ndim
4
PyArray_Descr * descr ; /* Pointer to type struct */
Common types in include int
, float
, bool
:
np.array ([-1, 0, 1], dtype=int)
array([-1, 0, 1])
np.array([-1, 0, 1], dtype=float)
array([-1., 0., 1.])
np.array([-1 , 0, 1], dtype=bool)
array([ True, False, True], 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)
array([(0, True), (1, False)], dtype=[('value', '<i8'), ('status', '?')])
This is called a structured array.
npy_intp *strides; /* bytes to jump to get */
/* to the next element */
x = np.arange(12).reshape((3, 4))
x
array([[ 0, 1, 2, 3], [ 4, 5, 6, 7], [ 8, 9, 10, 11]])
x.dtype
dtype('int64')
x.dtype.itemsize
8
x.strides
(32, 8)
We see (4 * itemsize, itemsize)
or (skip_bytes_row, skip_bytes_col)
.
We can take the transpose of the array without copying any data by simply manipulating shape
and strides
. What should they be?
int flags; /* Flags */
x = np.array([1, 2, 3])
z = x[::2]
x.flags
C_CONTIGUOUS : True F_CONTIGUOUS : True OWNDATA : True WRITEABLE : True ALIGNED : True UPDATEIFCOPY : False
z.flags
C_CONTIGUOUS : False F_CONTIGUOUS : False OWNDATA : False WRITEABLE : True ALIGNED : True UPDATEIFCOPY : False
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)
array([(0, True), (1, False)], dtype=[('value', '<i8'), ('status', '?')])
This is called a structured array, and is accessed like a dictionary:
x = np.array([(0, True), (1, False)], dtype=dt)
x['value']
array([0, 1])
x['status']
array([ True, False], dtype=bool)
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']
Combining of differently shaped arrays without creating large intermediate arrays:
x = np.arange(4)
print(x)
print(x + 3)
[0 1 2 3] [3 4 5 6]
See docstring of np.doc.broadcasting
for more detail.
a = np.arange(12).reshape((3, 4))
b = np.array([1, 2, 3])[:, np.newaxis]
print(a)
print()
print(b)
[[ 0 1 2 3] [ 4 5 6 7] [ 8 9 10 11]] [[1] [2] [3]]
a + b
array([[ 1, 2, 3, 4], [ 6, 7, 8, 9], [11, 12, 13, 14]])
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)
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) XXX
x = np.zeros((3, 5, 1))
y = np.zeros((3, 5, 8))
xx, yy = np.broadcast_arrays(x, y)
xx.shape
(3, 5, 8)
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
(3, 5, 8)
An ndarray can be indexed in two ways:
x = np.arange(9).reshape((3, 3))
x
array([[0, 1, 2], [3, 4, 5], [6, 7, 8]])
x[:, [1, 1, 2]]
array([[1, 1, 2], [4, 4, 5], [7, 7, 8]])
Which is equivalent to:
np.array((x[:, 1], x[:, 1], x[:, 2])).T
array([[1, 1, 2], [4, 4, 5], [7, 7, 8]])
Start with the following input:
x = np.array([[0, 1, 2], [3, 4, 5], [6, 7, 8]])
x.shape
(3, 3)
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]
idx0.shape, idx1.shape
((2, 2), (1, 2))
a, b = np.broadcast_arrays(idx0, idx1)
a
array([[0, 1], [1, 2]])
b
array([[0, 1], [0, 1]])
x
array([[0, 1, 2], [3, 4, 5], [6, 7, 8]])
x[idx0, idx1]
array([[0, 4], [3, 7]])
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)
(3, 2) (1, 2)
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!
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)
(3, 2) (1, 2)
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)
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 Question
Greetings,
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
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
array([[10, 5, 8, 8, 12, 14, 8, 10, 7, 9, 12, 11, 13, 13, 6], [11, 12, 12, 13, 6, 10, 14, 13, 14, 9, 8, 5, 8, 10, 5], [ 7, 8, 13, 6, 8, 8, 8, 12, 5, 6, 14, 14, 5, 9, 12], [ 8, 7, 12, 7, 5, 5, 9, 10, 10, 11, 13, 9, 6, 9, 14], [13, 6, 6, 12, 14, 14, 8, 11, 12, 7, 5, 8, 10, 14, 9], [ 9, 11, 9, 9, 8, 9, 9, 13, 9, 8, 12, 10, 10, 5, 6], [10, 14, 8, 5, 10, 5, 6, 7, 9, 7, 5, 8, 7, 5, 12], [10, 14, 5, 7, 12, 7, 14, 7, 8, 8, 7, 8, 9, 6, 7], [14, 6, 9, 11, 13, 7, 8, 5, 5, 11, 5, 11, 8, 8, 13], [13, 13, 7, 8, 7, 5, 13, 13, 8, 13, 7, 13, 9, 8, 5]])
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!
slices = block[idx_i, idx_j, idx_k]
slices.shape
(10, 15, 7)
(10, 15, 7)
Now verify that our window is centered on k
everywhere:
slices[:, :, 3]
array([[10, 5, 8, 8, 12, 14, 8, 10, 7, 9, 12, 11, 13, 13, 6], [11, 12, 12, 13, 6, 10, 14, 13, 14, 9, 8, 5, 8, 10, 5], [ 7, 8, 13, 6, 8, 8, 8, 12, 5, 6, 14, 14, 5, 9, 12], [ 8, 7, 12, 7, 5, 5, 9, 10, 10, 11, 13, 9, 6, 9, 14], [13, 6, 6, 12, 14, 14, 8, 11, 12, 7, 5, 8, 10, 14, 9], [ 9, 11, 9, 9, 8, 9, 9, 13, 9, 8, 12, 10, 10, 5, 6], [10, 14, 8, 5, 10, 5, 6, 7, 9, 7, 5, 8, 7, 5, 12], [10, 14, 5, 7, 12, 7, 14, 7, 8, 8, 7, 8, 9, 6, 7], [14, 6, 9, 11, 13, 7, 8, 5, 5, 11, 5, 11, 8, 8, 13], [13, 13, 7, 8, 7, 5, 13, 13, 8, 13, 7, 13, 9, 8, 5]])
np.all(slices[:, :, 3] == k)
True
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):
data = np.array([1, 2, 3, 4])
print(data.shape, data.strides)
(4,) (8,)
N = 400000
repeat = np.lib.stride_tricks.as_strided(data, shape=(N, 4), strides=(0, 8))
repeat
array([[1, 2, 3, 4], [1, 2, 3, 4], [1, 2, 3, 4], ..., [1, 2, 3, 4], [1, 2, 3, 4], [1, 2, 3, 4]])
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, 1, 9)
(2 * 28) + (1 * 30) + (9 * 31)
365
solve_diophantine(month_days, (1, 20, 20), 365)
(1, 4, 7)
(1 * 28) + (4 * 30) + (7 * 31)
365
%reload_ext load_style
%load_style css/notebook.css