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

https://github.com/rougier/numpy-tutorial

Setup

In [2]:
import numpy as np  # We always use this convention,
                    # also in the documentation
In [3]:
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.

In [5]:
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

The NumPy ndarray

Revision: Structure of an array

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:

In [6]:
x = np.array([1, 2, 3])
In [7]:
x.dtype
Out[7]:
dtype('int64')
In [10]:
x.ctypes.data
Out[10]:
18269632
In [14]:
import ctypes
list(ctypes.string_at(x.ctypes.data, x.size * x.itemsize))
Out[14]:
[1, 0, 0, 0, 0, 0, 0, 0, 2, 0, 0, 0, 0, 0, 0, 0, 3, 0, 0, 0, 0, 0, 0, 0]

Dimensions

int nd;
npy_intp *dimensions; /* number of dimensions */
                      /* size in each dimension */


In [15]:
x = np.array([])
In [16]:
np.array(0).shape
Out[16]:
()
In [17]:
x = np.random.random((3, 2, 3, 3))
x.shape
Out[17]:
(3, 2, 3, 3)
In [18]:
x.ndim
Out[18]:
4

Data type descriptors

PyArray_Descr * descr ;  /* Pointer to type struct */

Common types in include int, float, bool:

In [19]:
np.array ([-1, 0, 1], dtype=int)
Out[19]:
array([-1,  0,  1])
In [20]:
np.array([-1, 0, 1], dtype=float)
Out[20]:
array([-1.,  0.,  1.])
In [21]:
np.array([-1 , 0, 1], dtype=bool)
Out[21]:
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:

In [22]:
dt = np.dtype([('value', np.int), ('status', np.bool)])
np.array([(0, True), (1, False)], dtype=dt)
Out[22]:
array([(0, True), (1, False)], 
      dtype=[('value', '<i8'), ('status', '?')])

This is called a structured array.

Strides

npy_intp *strides;  /* bytes to jump to get */
                    /* to the next element */
In [23]:
x = np.arange(12).reshape((3, 4))
x
Out[23]:
array([[ 0,  1,  2,  3],
       [ 4,  5,  6,  7],
       [ 8,  9, 10, 11]])
In [24]:
x.dtype
Out[24]:
dtype('int64')
In [25]:
x.dtype.itemsize
Out[25]:
8
In [26]:
x.strides
Out[26]:
(32, 8)

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 */
In [27]:
x = np.array([1, 2, 3])
z = x[::2]

x.flags
Out[27]:
  C_CONTIGUOUS : True
  F_CONTIGUOUS : True
  OWNDATA : True
  WRITEABLE : True
  ALIGNED : True
  UPDATEIFCOPY : False
In [28]:
z.flags
Out[28]:
  C_CONTIGUOUS : False
  F_CONTIGUOUS : False
  OWNDATA : False
  WRITEABLE : True
  ALIGNED : True
  UPDATEIFCOPY : False

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:

In [29]:
dt = np.dtype([('value', np.int), ('status', np.bool)])
np.array([(0, True), (1, False)], dtype=dt)
Out[29]:
array([(0, True), (1, False)], 
      dtype=[('value', '<i8'), ('status', '?')])

This is called a structured array, and is accessed like a dictionary:

In [30]:
x = np.array([(0, True), (1, False)], dtype=dt)
x['value']
Out[30]:
array([0, 1])
In [31]:
x['status']
Out[31]:
array([ True, False], dtype=bool)

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:

In [33]:
x = np.arange(4)
print(x)
print(x + 3)
[0 1 2 3]
[3 4 5 6]

Broadcasting in 1D

See docstring of np.doc.broadcasting for more detail.

Broadcasting (2D)

In [36]:
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]]
In [37]:
a + b
Out[37]:
array([[ 1,  2,  3,  4],
       [ 6,  7,  8,  9],
       [11, 12, 13, 14]])

Broadcasting (3D)

In [38]:
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)       XXX

Explicit broadcasting

In [81]:
x = np.zeros((3, 5, 1))
y = np.zeros((3, 5, 8))

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

xx.shape
Out[81]:
(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]]
In [85]:
np.broadcast_to(x, (3,5,8)).shape
Out[85]:
(3, 5, 8)

Fancy indexing

An ndarray can be indexed in two ways:

  • Using slices and scalars
  • Using other ndarrays ("fancy indexing")
In [40]:
x = np.arange(9).reshape((3, 3))
x
Out[40]:
array([[0, 1, 2],
       [3, 4, 5],
       [6, 7, 8]])
In [41]:
x[:, [1, 1, 2]]
Out[41]:
array([[1, 1, 2],
       [4, 4, 5],
       [7, 7, 8]])

Which is equivalent to:

In [42]:
np.array((x[:, 1], x[:, 1], x[:, 2])).T
Out[42]:
array([[1, 1, 2],
       [4, 4, 5],
       [7, 7, 8]])

Output shape of an indexing op

  1. Broadcast all index arrays against one another.
  2. Use the dimensions of slices as-is.

Start with the following input:

In [43]:
x = np.array([[0, 1, 2], [3, 4, 5], [6, 7, 8]])
In [44]:
x.shape
Out[44]:
(3, 3)

What would happen if we do

In [45]:
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)

In [46]:
idx0.shape, idx1.shape
Out[46]:
((2, 2), (1, 2))
In [47]:
a, b = np.broadcast_arrays(idx0, idx1)
a
Out[47]:
array([[0, 1],
       [1, 2]])
In [48]:
b
Out[48]:
array([[0, 1],
       [0, 1]])
In [49]:
x
Out[49]:
array([[0, 1, 2],
       [3, 4, 5],
       [6, 7, 8]])
In [50]:
x[idx0, idx1]
Out[50]:
array([[0, 4],
       [3, 7]])

Output shape of an indexing op (cont'd)

In [51]:
x = np.random.random((15, 12, 16, 3))

index_one = np.array([[0, 1], [2, 3], [4, 5]])
index_two = np.array([[0, 1]])
In [53]:
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!

Output shape of an indexing op (cont'd)

In [54]:
x = np.random.random((15, 12, 16, 3))
index_one = np.array([[0, 1], [2, 3], [4, 5]])
index_two = np.array([[0, 1]])
In [56]:
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)

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 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

Jack's cube

Test setup

In [87]:
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.

In [58]:
block = np.empty((ni, nj, nk) , dtype=int)
block[:] = np.arange(nk)[np.newaxis, np.newaxis, :]

Pick out a random fake horizon in k:

In [88]:
rng = np.random.RandomState(0)

k = rng.randint(5, 15, size=(ni, nj))
k
Out[88]:
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]])

Solution

These two indices assure that we take a slice at each (i, j) position

In [67]:
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:

In [68]:
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

In [70]:
slices = block[idx_i, idx_j, idx_k]
slices.shape
(10, 15, 7)
Out[70]:
(10, 15, 7)

Now verify that our window is centered on k everywhere:

In [71]:
slices[:, :, 3]
Out[71]:
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]])
In [72]:
np.all(slices[:, :, 3] == k)
Out[72]:
True

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

In [74]:
data = np.array([1, 2, 3, 4])
print(data.shape, data.strides)
(4,) (8,)
In [75]:
N = 400000
repeat = np.lib.stride_tricks.as_strided(data, shape=(N, 4), strides=(0, 8))
In [76]:
repeat
Out[76]:
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]])

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?

In [102]:
from numpy.core.multiarray_tests import solve_diophantine

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

solve_diophantine(month_days, limits, 365)
Out[102]:
(2, 1, 9)
In [103]:
(2 * 28) + (1 * 30) + (9 * 31) 
Out[103]:
365
In [105]:
solve_diophantine(month_days, (1, 20, 20), 365)
Out[105]:
(1, 4, 7)
In [106]:
(1 * 28) + (4 * 30) + (7 * 31)
Out[106]:
365

Questions, discussion and exercises


In [ ]:
%reload_ext load_style
%load_style css/notebook.css