2.2. 高级NumPy

作者Pauli Virtanen

NumPy是Python的科学工具栈的基础。它的目的是对一块内存中的许多元素实现高效的操作。理解它如何工作的细节有助于有效地利用其灵活性、采用有用的快捷操作。

本节包括:

  • NumPy数组的剖析及其结果。提示和技巧。
  • 通用函数:什么是通用函数、为什么使用通用函数、以及如果你想要一个新的通用函数要做什么。
  • 与其他工具集成:NumPy提供了几种方法来在ndarray中包装任何数据,没有不必要的拷贝。
  • 最近添加的功能,以及其中的内容:PEP 3118缓冲区,广义ufuncs ...

先决条件

  • NumPy
  • Cython
  • Pillow(Python图像库,在几个例子中用到)

小费

在本节中,将导入numpy,如下所示:

>>> import numpy as np

2.2.1. ndarray的生命

2.2.1.1. 它是...

ndarray =

内存块 + 索引方案 + 数据类型描述符

  • 原始数据
  • 如何定位元素
  • 如何解释元素
../../_images/threefundamental.png
typedef struct PyArrayObject {
PyObject_HEAD
/* Block of memory */
char *data;
/* Data type descriptor */
PyArray_Descr *descr;
/* Indexing scheme */
int nd;
npy_intp *dimensions;
npy_intp *strides;
/* Other stuff */
PyObject *base;
int flags;
PyObject *weakreflist;
} PyArrayObject;

2.2.1.2. 内存块

>>> x = np.array([1, 2, 3], dtype=np.int32)
>>> x.data
<... at ...>
>>> str(x.data)
'\x01\x00\x00\x00\x02\x00\x00\x00\x03\x00\x00\x00'

数据的内存地址:

>>> x.__array_interface__['data'][0] 
64803824

整个__array_interface__

>>> x.__array_interface__  
{'data': (35828928, False),
'descr': [('', '<i4')],
'shape': (4,),
'strides': None,
'typestr': '<i4',
'version': 3}

提醒:两个ndarrays可能共享同一个内存:

>>> x = np.array([1, 2, 3, 4])
>>> y = x[:-1]
>>> x[0] = 9
>>> y
array([9, 2, 3])

内存不需要由ndarray拥有:

>>> x = b'1234'      # The 'b' is for "bytes", necessary in Python 3

x是一个字符串(在Python 3中是一个字节),我们可以将其数据表示为int数组:

>>> y = np.frombuffer(x, dtype=np.int8)
>>> y.data
<... at ...>
>>> y.base is x
True
>>> y.flags
C_CONTIGUOUS : True
F_CONTIGUOUS : True
OWNDATA : False
WRITEABLE : False
ALIGNED : True
UPDATEIFCOPY : False

owndatawriteable标志指示内存块的状态。

2.2.1.3. 数据类型

2.2.1.3.1. 描述符

dtype描述数组中的单个项:

type

数据的标量类型,下面中的一个:

int8、int16、float64等。(固定尺寸)

str、unicode、void(大小可变)

itemize size
byteorder 字节顺序: 大端序 > / 小端序 < / 不可用|
fields 子类型,如果它是结构化数据类型
shape 如果是一个子数组,则它为数组的形状
>>> np.dtype(int).type      
<type 'numpy.int64'>
>>> np.dtype(int).itemsize
8
>>> np.dtype(int).byteorder
'='

2.2.1.3.2. 示例:读取.wav文件

.wav文件头:

chunk_id "RIFF"
chunk_size 4字节无符号小端整数
format "WAVE"
fmt_id “fmt
fmt_size 4字节无符号小端整数
audio_fmt 2字节无符号小端整数
num_channels 2字节无符号小端整数
sample_rate 4字节无符号小端整数
byte_rate 4字节无符号小端整数
block_align 2字节无符号小端整数
bits_per_sample 2字节无符号小端整数
data_id "data"
data_size 4字节无符号小端整数
  • 44字节的原始数据块(在文件的开头)
  • ...后跟实际声音数据的data_size字节。

.wav文件头作为NumPy结构化数据类型:

>>> wav_header_dtype = np.dtype([
... ("chunk_id", (bytes, 4)), # flexible-sized scalar type, item size 4
... ("chunk_size", "<u4"), # little-endian unsigned 32-bit integer
... ("format", "S4"), # 4-byte string
... ("fmt_id", "S4"),
... ("fmt_size", "<u4"),
... ("audio_fmt", "<u2"), #
... ("num_channels", "<u2"), # .. more of the same ...
... ("sample_rate", "<u4"), #
... ("byte_rate", "<u4"),
... ("block_align", "<u2"),
... ("bits_per_sample", "<u2"),
... ("data_id", ("S1", (2, 2))), # sub-array, just for fun!
... ("data_size", "u4"),
... #
... # the sound data itself cannot be represented here:
... # it does not have a fixed size
... ])

另见

wavreader.py

>>> wav_header_dtype['format']
dtype('S4')
>>> wav_header_dtype.fields
dict_proxy({'block_align': (dtype('uint16'), 32), 'format': (dtype('S4'), 8), 'data_id': (dtype(('S1', (2, 2))), 36), 'fmt_id': (dtype('S4'), 12), 'byte_rate': (dtype('uint32'), 28), 'chunk_id': (dtype('S4'), 0), 'num_channels': (dtype('uint16'), 22), 'sample_rate': (dtype('uint32'), 24), 'bits_per_sample': (dtype('uint16'), 34), 'chunk_size': (dtype('uint32'), 4), 'fmt_size': (dtype('uint32'), 16), 'data_size': (dtype('uint32'), 40), 'audio_fmt': (dtype('uint16'), 20)})
>>> wav_header_dtype.fields['format']
(dtype('S4'), 8)
  • 第一个元素是结构化数据中的子类型,对应于名称format
  • 第二个是它从元素开始的偏移量(以字节为单位)

练习

迷你练习,通过使用偏移量生成一个“稀疏”的dtype,并只有一些字段:

>>> wav_header_dtype = np.dtype(dict(
... names=['format', 'sample_rate', 'data_id'],
... offsets=[offset_1, offset_2, offset_3], # counted from start of structure in bytes
... formats=list of dtypes for each of the fields,
... ))

并用它来读取采样率,data_id(作为子数组)。

>>> f = open('data/test.wav', 'r')
>>> wav_header = np.fromfile(f, dtype=wav_header_dtype, count=1)
>>> f.close()
>>> print(wav_header)
[ ('RIFF', 17402L, 'WAVE', 'fmt ', 16L, 1, 1, 16000L, 32000L, 2, 16, [['d', 'a'], ['t', 'a']], 17366L)]
>>> wav_header['sample_rate']
array([16000], dtype=uint32)

让我们尝试访问子数组:

>>> wav_header['data_id']  
array([[['d', 'a'],
['t', 'a']]],
dtype='|S1')
>>> wav_header.shape
(1,)
>>> wav_header['data_id'].shape
(1, 2, 2)

当访问子数组时,维度被添加到末尾!

注意

有已经存在的模块,如wavfileaudiolab用于加载声音数据...

2.2.1.3.3. 转换和重新解释/视图

转换

  • 在赋值的时候
  • 在数组构造的时候
  • 在算术计算的时候
  • 等等
  • 以及手动转换:.astype(dtype)

数据重新解释

  • 手动去做:.view(dtype)
2.2.1.3.3.1. Casting
  • 算术计算中的转换,简而言之:

    • 只有类型(不是值!)的操作数事项
    • 选择能够代表两者的最大“安全”类型
    • 在某些情况下,标量可以“丢失”到数组
  • 通常拷贝数据时的转换:

    >>> x = np.array([1, 2, 3, 4], dtype=np.float)
    
    >>> x
    array([ 1., 2., 3., 4.])
    >>> y = x.astype(np.int8)
    >>> y
    array([1, 2, 3, 4], dtype=int8)
    >>> y + 1
    array([2, 3, 4, 5], dtype=int8)
    >>> y + 256
    array([257, 258, 259, 260], dtype=int16)
    >>> y + 256.0
    array([ 257., 258., 259., 260.])
    >>> y + np.array([256], dtype=np.int32)
    array([257, 258, 259, 260], dtype=int32)
  • 在setitem时的类型转换:在元素赋值时数组的dtype不更改:

    >>> y[:] = y + 1.5
    
    >>> y
    array([2, 3, 4, 5], dtype=int8)

注意

确切规则:请参阅numpy文档

2.2.1.3.3.2. 重新解释/视图
  • 内存中的数据块(4字节)

    0x01 || 0x02 || 0x03 || 0x04
    • 4个uint8,或者
    • 4个int8,或者
    • 2个int16,或者
    • 1个int32,或者
    • 1个float32,或者
    • ...

    如何从一个切换到另一个?

  1. 切换dtype:

    >>> x = np.array([1, 2, 3, 4], dtype=np.uint8)
    
    >>> x.dtype = "<i2"
    >>> x
    array([ 513, 1027], dtype=int16)
    >>> 0x0201, 0x0403
    (513, 1027)
0x01 0x02 || 0x03 0x04

注意

little-endian:在内存中,最低有效字节位于左边

  1. 创建新视图:

    >>> y = x.view("<i4")
    
    >>> y
    array([67305985], dtype=int32)
    >>> 0x04030201
    67305985
0x01 0x02 0x03 0x04

注意

  • .view()生成视图,不复制(或更改)内存块

  • 只改变dtype(并调整数组形状):

    >>> x[1] = 5
    
    >>> y
    array([328193], dtype=int32)
    >>> y.base is x
    True

小练习:数据重新解释

另见

view-colors.py

数组中有RGBA数据:

>>> x = np.zeros((10, 10, 4), dtype=np.int8)
>>> x[:, :, 0] = 1
>>> x[:, :, 1] = 2
>>> x[:, :, 2] = 3
>>> x[:, :, 3] = 4

其中最后三个维度是R,B和G,以及α通道。

如何做一个字段名称为'r','g','b','a'的(10, 10)结构化数组而不复制数据?

>>> y = ...                     
>>> assert (y['r'] == 1).all()
>>> assert (y['g'] == 2).all()
>>> assert (y['b'] == 3).all()
>>> assert (y['a'] == 4).all()

答案

...

警告

另一个正好占用4个字节内存的数组:

>>> y = np.array([[1, 3], [2, 4]], dtype=np.uint8).transpose()
>>> x = y.copy()
>>> x
array([[1, 2],
[3, 4]], dtype=uint8)
>>> y
array([[1, 2],
[3, 4]], dtype=uint8)
>>> x.view(np.int16)
array([[ 513],
[1027]], dtype=int16)
>>> 0x0201, 0x0403
(513, 1027)
>>> y.view(np.int16)
array([[ 769, 1026]], dtype=int16)
  • 发生了什么?
  • ...我们需要看看x[0,1]实际意味着什么
>>> 0x0301, 0x0402
(769, 1026)

2.2.1.4. 索引方案:strides

2.2.1.4.1. 要点

问题

>>> x = np.array([[1, 2, 3],
... [4, 5, 6],
... [7, 8, 9]], dtype=np.int8)
>>> str(x.data)
'\x01\x02\x03\x04\x05\x06\x07\x08\t'
At which byte in ``x.data`` does the item ``x[1, 2]`` begin?

答案(在NumPy中)

  • strides:找到下一个元素跳转的字节数
  • 每个维度一个stride
  >>> x.strides
(3, 1)
>>> byte_offset = 3*1 + 1*2 # to find x[1, 2]
>>> x.flat[byte_offset]
6
>>> x[1, 2]
6
- simple, **flexible**
2.2.1.4.1.1. C和Fortran的顺序
>>> x = np.array([[1, 2, 3],
... [4, 5, 6]], dtype=np.int16, order='C')
>>> x.strides
(6, 2)
>>> str(x.data)
'\x01\x00\x02\x00\x03\x00\x04\x00\x05\x00\x06\x00'
  • 需要跳转6个字节才能找到下一行
  • 需要跳转2个字节才能找到下一列
>>> y = np.array(x, order='F')
>>> y.strides
(2, 4)
>>> str(y.data)
'\x01\x00\x04\x00\x02\x00\x05\x00\x03\x00\x06\x00'
  • 需要跳转2个字节才能找到下一行
  • 需要跳转4个字节才能找到下一个列
  • 更高的维度类似:

    • C:最后的维度变化最快(=较小的stride)
    • F:第一个维度变化最快

    \mathrm{shape} &= (d_1, d_2, ..., d_n)
\\
\mathrm{strides} &= (s_1, s_2, ..., s_n)
\\
s_j^C &= d_{j+1} d_{j+2} ... d_{n} \times \mathrm{itemsize}
\\
s_j^F &= d_{1} d_{2} ... d_{j-1} \times \mathrm{itemsize}

注意

现在我们可以理解.view()的行为:

>>> y = np.array([[1, 3], [2, 4]], dtype=np.uint8).transpose()
>>> x = y.copy()

转置不影响数据的内存布局,只影响stride

>>> x.strides
(2, 1)
>>> y.strides
(1, 2)
>>> str(x.data)  
'\x01\x02\x03\x04'
>>> str(y.data)
'\x01\x03\x02\x04'
  • 当解释为2个int16时,结果不同
  • .copy()以C顺序创建新数组(默认情况下)
2.2.1.4.1.2. 用整数切片
  • 所有东西都可以通过只更改shapestrides来表示,并可能只通过调整data指针表示!
  • 永远不会复制数据
>>> x = np.array([1, 2, 3, 4, 5, 6], dtype=np.int32)
>>> y = x[::-1]
>>> y
array([6, 5, 4, 3, 2, 1], dtype=int32)
>>> y.strides
(-4,)
>>> y = x[2:]
>>> y.__array_interface__['data'][0] - x.__array_interface__['data'][0]
8
>>> x = np.zeros((10, 10, 10), dtype=np.float)
>>> x.strides
(800, 80, 8)
>>> x[::2,::3,::4].strides
(1600, 240, 32)
  • 类似地,转置永远不会生成拷贝(它只是交换stride):

    >>> x = np.zeros((10, 10, 10), dtype=np.float)
    
    >>> x.strides
    (800, 80, 8)
    >>> x.T.strides
    (8, 80, 800)

但是:不是所有的reshape操作都可以通过使用stride来表示:

>>> a = np.arange(6, dtype=np.int8).reshape(3, 2)
>>> b = a.T
>>> b.strides
(1, 2)

到现在为止还挺好。然而:

>>> str(a.data)  
'\x00\x01\x02\x03\x04\x05'
>>> b
array([[0, 2, 4],
[1, 3, 5]], dtype=int8)
>>> c = b.reshape(3*2)
>>> c
array([0, 2, 4, 1, 3, 5], dtype=int8)

这里,只用一个stride和a的内存块,没有办法表示数组c因此,reshape操作需要在此处进行复制。

2.2.1.4.2. Example: fake dimensions with strides

Stride操纵

>>> from numpy.lib.stride_tricks import as_strided
>>> help(as_strided)
as_strided(x, shape=None, strides=None)
Make an ndarray from the given array with the given shape and strides

警告

as_strided不会检查是否在内存块边界内...

>>> x = np.array([1, 2, 3, 4], dtype=np.int16)
>>> as_strided(x, strides=(2*2, ), shape=(2, ))
array([1, 3], dtype=int16)
>>> x[::2]
array([1, 3], dtype=int16)

另见

stride-fakedims.py

练习

array([1, 2, 3, 4], dtype=np.int8)
-> array([[1, 2, 3, 4],
[1, 2, 3, 4],
[1, 2, 3, 4]], dtype=np.int8)

只使用as_strided

Hint: byte_offset = stride[0]*index[0] + stride[1]*index[1] + ...

Spoiler

...

2.2.1.4.3. Broadcasting

  • 用它做一些有用的东西:[1, 2, 3, 4][5, 6, 7]的外积
>>> x = np.array([1, 2, 3, 4], dtype=np.int16)
>>> x2 = as_strided(x, strides=(0, 1*2), shape=(3, 4))
>>> x2
array([[1, 2, 3, 4],
[1, 2, 3, 4],
[1, 2, 3, 4]], dtype=int16)
>>> y = np.array([5, 6, 7], dtype=np.int16)
>>> y2 = as_strided(y, strides=(1*2, 0), shape=(3, 4))
>>> y2
array([[5, 5, 5, 5],
[6, 6, 6, 6],
[7, 7, 7, 7]], dtype=int16)
>>> x2 * y2
array([[ 5, 10, 15, 20],
[ 6, 12, 18, 24],
[ 7, 14, 21, 28]], dtype=int16)

...似乎有点熟悉...

>>> x = np.array([1, 2, 3, 4], dtype=np.int16)
>>> y = np.array([5, 6, 7], dtype=np.int16)
>>> x[np.newaxis,:] * y[:,np.newaxis]
array([[ 5, 10, 15, 20],
[ 6, 12, 18, 24],
[ 7, 14, 21, 28]], dtype=int16)
  • 在内部,数组广播的确使用0步幅实现。

2.2.1.4.4. More tricks: diagonals

另见

stride-diagonals.py

挑战

  • 选择矩阵的对角项:(假设C存储顺序):

    >>> x = np.array([[1, 2, 3],
    
    ... [4, 5, 6],
    ... [7, 8, 9]], dtype=np.int32)
    >>> x_diag = as_strided(x, shape=(3,), strides=(???,))
  • 选择第一个超对角线条目[2, 6]

  • 和子对角线?

(提示到最后两个:切片首先移动点大步
从。开始)。

答案

...

另见

stride-diagonals.py

挑战

计算张量图:

>>> x = np.arange(5*5*5*5).reshape(5, 5, 5, 5)
>>> s = 0
>>> for i in range(5):
... for j in range(5):
... s += x[j, i, j, i]

,并对结果使用sum()

>>> y = as_strided(x, shape=(5, 5), strides=(TODO, TODO))   
>>> s2 = ...
>>> assert s == s2

答案

...

2.2.1.4.5. CPU缓存的影响

内存布局可能会影响性能:

In [1]: x = np.zeros((20000,))
In [2]: y = np.zeros((20000*67,))[::67]
In [3]: x.shape, y.shape
((20000,), (20000,))
In [4]: %timeit x.sum()
100000 loops, best of 3: 0.180 ms per loop
In [5]: %timeit y.sum()
100000 loops, best of 3: 2.34 ms per loop
In [6]: x.strides, y.strides
((8,), (536,))

较小的步幅更快?

../../_images/cpu-cacheline.png
  • CPU以块为单位将数据从主存储器拉到其缓存
  • 如果连续操作的许多数组元素在一个单一的块内(stride较小):
    • \Rightarrow减少转移所需
    • \Rightarrow更快

另见

numexpr旨在减轻数组计算中的缓存影响。

2.2.1.4.6. Example: inplace operations (caveat emptor)

  • 有时,

    >>> a -= b    
    

    与下面的语句是不等同的

    >>> a -= b.copy()    
    
>>> x = np.array([[1, 2], [3, 4]])
>>> x -= x.transpose()
>>> x
array([[ 0, -1],
[ 4, 0]])
>>> y = np.array([[1, 2], [3, 4]])
>>> y -= y.T.copy()
>>> y
array([[ 0, -1],
[ 1, 0]])
  • xx.transpose()共享数据
  • x -= x.transpose()逐个元素修改数据...
  • 因为xx.transpose()有不同的striding,修改的数据重新出现在RHS

2.2.1.5. 剖析结果

../../_images/threefundamental.png
  • 内存块:可以共享,.base.data
  • data type descriptor: structured data, sub-arrays, byte order, casting, viewing, .astype(), .view()
  • strided indexing: strides, C/F-order, slicing w/ integers, as_strided, broadcasting, stride tricks, diag, CPU cache coherence

2.2.2. 通用函数

2.2.2.1. 它们是什么?

  • Ufunc对数组的所有元素执行元素级别操作。

    例子:

    np.add, np.subtract, scipy.special.*, ...
    
  • 自动支持:broadcasting、casting...

  • ufunc的作者只需要提供元素级别操作,NumPy负责其余的操作。

  • 元素级别操作需要在C(或例如Cython)中实现,

2.2.2.1.1. Parts of an Ufunc

  1. 由用户提供

    void ufunc_loop(void **args, int *dimensions, int *steps, void *data)
    
    {
    /*
    * int8 output = elementwise_function(int8 input_1, int8 input_2)
    *
    * This function must compute the ufunc for many values at once,
    * in the way shown below.
    */
    char *input_1 = (char*)args[0];
    char *input_2 = (char*)args[1];
    char *output = (char*)args[2];
    int i;
    for (i = 0; i < dimensions[0]; ++i) {
    *output = elementwise_function(*input_1, *input_2);
    input_1 += steps[0];
    input_2 += steps[1];
    output += steps[2];
    }
    }
  2. NumPy部分,由

    char types[3]
    
    types[0] = NPY_BYTE /* type of first input arg */
    types[1] = NPY_BYTE /* type of second input arg */
    types[2] = NPY_BYTE /* type of third input arg */
    PyObject *python_ufunc = PyUFunc_FromFuncAndData(
    ufunc_loop,
    NULL,
    types,
    1, /* ntypes */
    2, /* num_inputs */
    1, /* num_outputs */
    identity_element,
    name,
    docstring,
    unused)
    • ufunc也可以支持多种不同输入—输出类型组合。

2.2.2.1.2. Making it easier

  1. ufunc_loop是非常通用的形式,NumPy提供预制的

    PyUfunc_f_f float elementwise_func(float input_1)
    PyUfunc_ff_f float elementwise_func(float input_1, float input_2)
    PyUfunc_d_d double elementwise_func(double input_1)
    PyUfunc_dd_d double elementwise_func(double input_1, double input_2) t0 >
    PyUfunc_D_D elementwise_func(npy_cdouble *input, npy_cdouble* output)
    PyUfunc_DD_D elementwise_func(npy_cdouble * in1, npy_cdouble * in2, npy_cdouble * out)
    • 只需提供elementwise_func
    • ...除非你的元素级别函数不在上述形式之一

2.2.2.2. 练习:从头开始构建ufunc

Mandelbrot分形由迭代定义

z \leftarrow z^2 + c

其中c = x + i y是一个复数。这个迭代被重复 - 如果z无论迭代运行多长时间,c属于Mandelbrot集。

  • 使ufunc称为mandel(z0, c)

    z = z0
    
    for k in range(iterations):
    z = z*z + c

    例如迭代100次,或者直到z.real**2 + z.imag**2 > 1000使用它来确定哪个c在Mandelbrot集合中。

  • 我们的函数是一个简单的函数,所以使用PyUFunc_* helpers。

  • 用Cython编写

另见

mandel.py,mandelplot.py

提醒:一些预制的Ufunc循环:

PyUfunc_f_f float elementwise_func(float input_1)
PyUfunc_ff_f float elementwise_func(float input_1, 浮点数 input_2) t0>
PyUfunc_d_d double elementwise_func(double input_1)
PyUfunc_dd_d double elementwise_func(double input_1, double input_2) t0 >
PyUfunc_D_D elementwise_func(complex_double * input, complex_double * 输出)
PyUfunc_DD_D elementwise_func(complex_double * in1, complex_double * in2, complex_double * out)

类型代码:

NPY_BOOL, NPY_BYTE, NPY_UBYTE, NPY_SHORT, NPY_USHORT, NPY_INT, NPY_UINT,
NPY_LONG, NPY_ULONG, NPY_LONGLONG, NPY_ULONGLONG, NPY_FLOAT, NPY_DOUBLE,
NPY_LONGDOUBLE, NPY_CFLOAT, NPY_CDOUBLE, NPY_CLONGDOUBLE, NPY_DATETIME,
NPY_TIMEDELTA, NPY_OBJECT, NPY_STRING, NPY_UNICODE, NPY_VOID

2.2.2.3. 解决方案:从头构建ufunc

# The elementwise function
# ------------------------
cdef void mandel_single_point(double complex *z_in,
double complex *c_in,
double complex *z_out) nogil:
#
# The Mandelbrot iteration
#
#
# Some points of note:
#
# - It's *NOT* allowed to call any Python functions here.
#
# The Ufunc loop runs with the Python Global Interpreter Lock released.
# Hence, the ``nogil``.
#
# - And so all local variables must be declared with ``cdef``
#
# - Note also that this function receives *pointers* to the data;
# the "traditional" solution to passing complex variables around
#
cdef double complex z = z_in[0]
cdef double complex c = c_in[0]
cdef int k # the integer we use in the for loop
# Straightforward iteration
for k in range(100):
z = z*z + c
if z.real**2 + z.imag**2 > 1000:
break
# Return the answer for this point
z_out[0] = z
# Boilerplate Cython definitions
#
# Pulls definitions from the Numpy C headers.
# -------------------------------------------
from numpy cimport import_array, import_ufunc
from numpy cimport (PyUFunc_FromFuncAndData,
PyUFuncGenericFunction)
from numpy cimport NPY_CDOUBLE
from numpy cimport PyUFunc_DD_D
# Required module initialization
# ------------------------------
import_array()
import_ufunc()
# The actual ufunc declaration
# ----------------------------
cdef PyUFuncGenericFunction loop_func[1]
cdef char input_output_types[3]
cdef void *elementwise_funcs[1]
loop_func[0] = PyUFunc_DD_D
input_output_types[0] = NPY_CDOUBLE
input_output_types[1] = NPY_CDOUBLE
input_output_types[2] = NPY_CDOUBLE
elementwise_funcs[0] = <void*>mandel_single_point
mandel = PyUFunc_FromFuncAndData(
loop_func,
elementwise_funcs,
input_output_types,
1, # number of supported input types
2, # number of input args
1, # number of output args
0, # `identity` element, never mind this
"mandel", # function name
"mandel(z, c) -> computes iterated z*z + c", # docstring
0 # unused
)
"""
Plot Mandelbrot
================
Plot the Mandelbrot ensemble.
"""
import numpy as np
import mandel
x = np.linspace(-1.7, 0.6, 1000)
y = np.linspace(-1.4, 1.4, 1000)
c = x[None,:] + 1j*y[:,None]
z = mandel.mandel(c, c)
import matplotlib.pyplot as plt
plt.imshow(abs(z)**2 < 1000, extent=[-1.7, 0.6, -1.4, 1.4])
plt.gray()
plt.show()
../../_images/mandelbrot.png

注意

大多数样板可以通过这些Cython模块自动化:

http://wiki.cython.org/MarkLodato/CreatingUfuncs

几个接受的输入类型

例如,支持单精度和双精度版本

cdef void mandel_single_point(double complex *z_in,
double complex *c_in,
double complex *z_out) nogil:
...
cdef void mandel_single_point_singleprec(float complex *z_in,
float complex *c_in,
float complex *z_out) nogil:
...
cdef PyUFuncGenericFunction loop_funcs[2]
cdef char input_output_types[3*2]
cdef void *elementwise_funcs[1*2]
loop_funcs[0] = PyUFunc_DD_D
input_output_types[0] = NPY_CDOUBLE
input_output_types[1] = NPY_CDOUBLE
input_output_types[2] = NPY_CDOUBLE
elementwise_funcs[0] = <void*>mandel_single_point
loop_funcs[1] = PyUFunc_FF_F
input_output_types[3] = NPY_CFLOAT
input_output_types[4] = NPY_CFLOAT
input_output_types[5] = NPY_CFLOAT
elementwise_funcs[1] = <void*>mandel_single_point_singleprec
mandel = PyUFunc_FromFuncAndData(
loop_func,
elementwise_funcs,
input_output_types,
2, # number of supported input types <----------------
2, # number of input args
1, # number of output args
0, # `identity` element, never mind this
"mandel", # function name
"mandel(z, c) -> computes iterated z*z + c", # docstring
0 # unused
)

2.2.2.4. 广义ufuncs

ufunc

output = elementwise_function(input)

outputinput只能是单个数组元素。

通用ufunc

outputinput可以是具有固定维数的数组

例如,矩阵的积(对角元素的总和):

input shape = (n, n)
output shape = () i.e. scalar
(n, n) -> ()

矩阵的积:

input_1 shape = (m, n)
input_2 shape = (n, p)
output shape = (m, p)
(m, n), (n, p) -> (m, p)
  • 这被称为广义ufunc的“签名”
  • g-ufunc所作用的维度是“核心维度”

NumPy中的状态

  • g-ufuncs在NumPy已经...

  • 新的可以用PyUFunc_FromFuncAndDataAndSignature

  • 大多数线性代数函数被实现为g-ufuncs以支持使用堆栈阵列:

    >>> import numpy as np
    
    >>> np.linalg.det(np.random.rand(3, 5, 5))
    array([ 0.00965823, -0.13344729, 0.04583961])
    >>> np.linalg._umath_linalg.det.signature
    '(m,m)->()'
  • 我们还提供了几个g-ufuncs进行测试,ATM:

    >>> import numpy.core.umath_tests as ut
    
    >>> ut.matrix_multiply.signature
    '(m,n),(n,p)->(m,p)'
    >>> x = np.ones((10, 2, 4))
    >>> y = np.ones((10, 4, 5))
    >>> ut.matrix_multiply(x, y).shape
    (10, 2, 5)
  • 在这两个示例中,最后两个维度变为核心维度,并根据签名

  • 否则,g-ufunc运行“元素级别”

  • 这种方式的矩阵乘法可以用于一次对许多小矩阵进行操作

广义ufunc循环

矩阵乘法(m,n),(n,p) - (m,p)

void gufunc_loop(void **args, int *dimensions, int *steps, void *data)
{
char *input_1 = (char*)args[0]; /* these are as previously */
char *input_2 = (char*)args[1];
char *output = (char*)args[2];
int input_1_stride_m = steps[3]; /* strides for the core dimensions */
int input_1_stride_n = steps[4]; /* are added after the non-core */
int input_2_strides_n = steps[5]; /* steps */
int input_2_strides_p = steps[6];
int output_strides_n = steps[7];
int output_strides_p = steps[8];
int m = dimension[1]; /* core dimensions are added after */
int n = dimension[2]; /* the main dimension; order as in */
int p = dimension[3]; /* signature */
int i;
for (i = 0; i < dimensions[0]; ++i) {
matmul_for_strided_matrices(input_1, input_2, output,
strides for each array...);
input_1 += steps[0];
input_2 += steps[1];
output += steps[2];
}
}

2.2.3. 互操作性功能

2.2.3.1. 共享多维,类型化数据

假设你

  1. 写一个库比句柄(多维)二进制数据,
  2. 想要使它容易操作数据与NumPy,或任何其他库,
  3. ...但是不是喜欢将NumPy作为依赖。

目前,3个解决方案:

  1. “旧”缓冲接口
  2. 数组接口
  3. “新”缓冲接口( PEP 3118

2.2.3.2. 旧缓冲区协议

  • 只有一维缓冲器
  • 无数据类型信息
  • C级接口; PyBufferProcs tp_as_buffer
  • 但它被集成到Python(例如字符串支持它)

使用枕头(Python成像库)的微型锻炼:

也可以看看

pilbuffer.py

>>> from PIL import Image
>>> data = np.zeros((200, 200, 4), dtype=np.int8)
>>> data[:, :] = [255, 0, 0, 255] # Red
>>> # In PIL, RGBA images consist of 32-bit integers whose bytes are [RR,GG,BB,AA]
>>> data = data.view(np.int32).squeeze()
>>> img = Image.frombuffer("RGBA", (200, 200), data, "raw", "RGBA", 0, 1)
>>> img.save('test.png')

Q:

检查data现在修改后,如果再次保存img,会发生什么情况。

2.2.3.3. 旧缓冲区协议

"""
From buffer
============
Show how to exchange data between numpy and a library that only knows
the buffer interface.
"""
import numpy as np
import Image
# Let's make a sample image, RGBA format
x = np.zeros((200, 200, 4), dtype=np.int8)
x[:,:,0] = 254 # red
x[:,:,3] = 255 # opaque
data = x.view(np.int32) # Check that you understand why this is OK!
img = Image.frombuffer("RGBA", (200, 200), data)
img.save('test.png')
#
# Modify the original data, and save again.
#
# It turns out that PIL, which knows next to nothing about Numpy,
# happily shares the same data.
#
x[:,:,1] = 254
img.save('test2.png')
../../_images/test.png ../../_images/test2.png

2.2.3.4. 数组接口协议

  • 多维缓冲区
  • 数据类型信息存在
  • NumPy特定方法;缓慢地弃用(但不会消失)
  • 没有集成在Python中
>>> x = np.array([[1, 2], [3, 4]])
>>> x.__array_interface__
{'data': (171694552, False), # memory address of data, is readonly?
'descr': [('', '<i4')], # data type descriptor
'typestr': '<i4', # same, in another form
'strides': None, # strides; or None if in C-order
'shape': (2, 2),
'version': 3,
}
::
>>> from PIL import Image
>>> img = Image.open('data/test.png')
>>> img.__array_interface__
{'data': ...,
'shape': (200, 200, 4),
'typestr': '|u1'}
>>> x = np.asarray(img)
>>> x.shape
(200, 200, 4)

注意

还定义了一个更友好的数组接口变体。

2.2.4. Array siblings:chararraymaskedarraymatrix

2.2.4.1. chararray:向量化字符串操作

>>> x = np.array(['a', '  bbb', '  ccc']).view(np.chararray)
>>> x.lstrip(' ')
chararray(['a', 'bbb', 'ccc'],
dtype='...')
>>> x.upper()
chararray(['A', ' BBB', ' CCC'],
dtype='...')

注意

.view()有第二个含义:它可以使ndarray成为一个专门的ndarray子类的实例

2.2.4.2. masked_array缺少数据

掩码数组是可能缺少或无效条目的数组。

例如,假设我们有一个数组,其中第四个条目无效:

>>> x = np.array([1, 2, 3, -99, 5])

描述这一点的一种方法是创建一个屏蔽数组:

>>> mx = np.ma.masked_array(x, mask=[0, 0, 0, 1, 0])
>>> mx
masked_array(data = [1 2 3 -- 5],
mask = [False False False True False],
fill_value = 999999)

屏蔽平均值忽略屏蔽数据:

>>> mx.mean()
2.75
>>> np.mean(mx)
2.75

警告

不是所有NumPy函数都会考虑掩码,例如np.dot,因此请检查返回类型。

masked_array返回视图到原始数组:

>>> mx[1] = 9
>>> x
array([ 1, 9, 3, -99, 5])

2.2.4.2.1. The mask

您可以通过以下方式修改掩码:

>>> mx[1] = np.ma.masked
>>> mx
masked_array(data = [1 -- 3 -- 5],
mask = [False True False True False],
fill_value = 999999)

掩码在分配时被清除:

>>> mx[1] = 9
>>> mx
masked_array(data = [1 9 3 -- 5],
mask = [False False False True False],
fill_value = 999999)

面罩也可直接使用:

>>> mx.mask
array([False, False, False, True, False], dtype=bool)

被掩码的条目可以用给定值填充以得到通常的阵列:

>>> x2 = mx.filled(-1)
>>> x2
array([ 1, 9, 3, -1, 5])

也可以清除掩码:

>>> mx.mask = np.ma.nomask
>>> mx
masked_array(data = [1 9 3 -99 5],
mask = [False False False False False],
fill_value = 999999)

2.2.4.2.2. Domain-aware functions

掩码数组包还包含域感知函数:

>>> np.ma.log(np.array([1, 2, -1, -2, 3, -5]))
masked_array(data = [0.0 0.6931471805599453 -- -- 1.0986122886681098 --],
mask = [False False True True False True],
fill_value = 1e+20)

注意

简化和更无缝的支持处理数组中缺失的数据正在进入NumPy 1.7。敬请关注!

示例:掩码统计

加拿大护林员在1903-1910年和1917-1918年计算野兔和ly were时分心,得到的数字是错误的。(胡萝卜农民保持警觉,虽然。)计算随时间的平均群体,忽略无效数字。

>>> data = np.loadtxt('data/populations.txt')
>>> populations = np.ma.masked_array(data[:,1:])
>>> year = data[:, 0]
>>> bad_years = (((year >= 1903) & (year <= 1910))
... | ((year >= 1917) & (year <= 1918)))
>>> # '&' means 'and' and '|' means 'or'
>>> populations[bad_years, 0] = np.ma.masked
>>> populations[bad_years, 1] = np.ma.masked
>>> populations.mean(axis=0)
masked_array(data = [40472.72727272727 18627.272727272728 42400.0],
mask = [False False False],
fill_value = 1e+20)
>>> populations.std(axis=0)
masked_array(data = [21087.656489006717 15625.799814240254 3322.5062255844787],
mask = [False False False],
fill_value = 1e+20)

注意Matplotlib知道掩码数组:

>>> plt.plot(year, populations, 'o-')   
[<matplotlib.lines.Line2D object at ...>, ...]

[源代码hires.pngpdf]

../../_images/numpy_intro_8.png

2.2.4.3. recarray:纯粹方便

>>> arr = np.array([('a', 1), ('b', 2)], dtype=[('x', 'S1'), ('y', int)])
>>> arr2 = arr.view(np.recarray)
>>> arr2.x
chararray(['a', 'b'],
dtype='|S1')
>>> arr2.y
array([1, 2])

2.2.4.4. matrix:方便?

  • 总是2-D
  • *是矩阵乘积,而不是元素级别
>>> np.matrix([[1, 0], [0, 1]]) * np.matrix([[1, 2], [3, 4]])
matrix([[1, 2],
[3, 4]])

2.2.5. 摘要

  • ndarray的解剖:data,dtype,strides。
  • 通用功能:元素级操作,如何制作新的
  • Ndarray子类
  • 各种缓冲接口,用于与其他工具集成
  • 最近增加:PEP 3118,广义ufuncs

2.2.6. 贡献NumPy/Scipy

2.2.6.1. 为什么

  • “有一个错误?
  • “我不明白这是应该做什么?
  • “我有这个花哨的代码。你想要吗?
  • “我想帮助!我能做什么?”

2.2.6.2. 报告错误

2.2.6.2.1. Good bug report

Title: numpy.random.permutations fails for non-integer arguments
I'm trying to generate random permutations, using numpy.random.permutations
When calling numpy.random.permutation with non-integer arguments
it fails with a cryptic error message::
>>> np.random.permutation(12)
array([11, 5, 8, 4, 6, 1, 9, 3, 7, 2, 10, 0])
>>> np.random.permutation(12.)
Traceback (most recent call last):
File "<stdin>", line 1, in <module>
File "mtrand.pyx", line 3311, in mtrand.RandomState.permutation
File "mtrand.pyx", line 3254, in mtrand.RandomState.shuffle
TypeError: len() of unsized object
This also happens with long arguments, and so
np.random.permutation(X.shape[0]) where X is an array fails on 64
bit windows (where shape is a tuple of longs).
It would be great if it could cast to integer or at least raise a
proper error for non-integer types.
I'm using NumPy 1.4.1, built from the official tarball, on Windows
64 with Visual studio 2008, on Python.org 64-bit Python.
  1. 你想做什么?

  2. 重现错误的小代码段(如果可能)

    • 实际发生了什么
    • 你会期望什么
  3. 平台(Windows / Linux / OSX,32/64位,x86 / PPC,...)

  4. NumPy / Scipy的版本

    >>> print(np.__version__) 
    
    1...

    检查以下是您期望的

    >>> print(np.__file__) 
    
    /...

    如果你有旧的/破碎的NumPy安装躺在。

    如果不确定,请尝试删除现有的NumPy安装,然后重新安装...

2.2.6.3. 参与文档

  1. 文档编辑器

    • http://docs.scipy.org/doc/numpy

    • 注册

      • 注册一个帐号

      • 订阅scipy-dev邮寄名单(仅限订阅者)

      • 邮件列表的问题:你收到邮件

      • 发送邮件@ scipy-dev邮件列表;请求激活:

        To: scipy-dev@scipy.org
        
        Hi,
        I'd like to edit NumPy/Scipy docstrings. My account is XXXXX
        Cheers,
        N. N.
  2. 编辑来源并发送修补程序(对于错误)

  3. 投诉邮件列表

2.2.6.4. 贡献功能

  1. 在邮件列表上询问,如果不确定它应该去哪里

  2. 编写补丁,在错误跟踪集上添加增强票

  3. OR,创建一个实现功能+添加增强票的Git分支。

    # Clone numpy repository
    
    git clone --origin svn http://projects.scipy.org/git/numpy.git numpy
    cd numpy
    # Create a feature branch
    git checkout -b name-of-my-feature-branch svn/trunk
    <edit stuff>
    git commit -a
    • https://github.com(或任何地方)创建帐户
    • 创建一个新的仓库@ Github
    • 将您的工作推送到github
    git remote add github git@github:USERNAME/REPOSITORYNAME.git
    
    git push github name-of-my-feature-branch

2.2.6.5. 如何帮助,一般

  • 错误修复总是欢迎!

    • 你最害怕什么
    • 浏览跟踪器
  • 文档工作

  • 在通信频道上询问:

    • numpy-discussion列表
    • scipy-dev列表