本篇重点介绍广播机制以及针对高维数组的轴操作介绍

本文重点介绍高维数组的广播机制和轴运算，最后介绍NumPy的C语言扩展。

1 次广播

NumPy 操作通常在两个数组的元素级别执行。在最简单的情况下，两个形状完全相同的数组操作，如下例所示，

a = np.array([1.0, 2.0, 3.0])

b = np.array([2.0, 2.0, 2.0])

a * b

Numpy 的广播机制是指在执行算术运算时处理不同形状的数组的方式。在某些规则下，较小的数组会在较大的数组上进行广播，以便数组具有兼容的形状。

a = np.array([1.0, 2.0, 3.0])

b = 2.0

a * b

发现这两个计算的结果是一样的，但是第二个是有广播机制在起作用。广播规则

当对两个数组执行操作时，NumPy 会比较它们的形状。它从形状的最右侧开始向左侧进行逐一比较。如果在以下两种情况之一中比较所有席位，

然后就可以对这两个数组进行操作了。如果不满足这些条件，则会引发 ValueError，指示数组具有不兼容的形状。

可以看出，数组的形状就像一个人的性格。如果两个人的性格不匹配，就不能在一起。

在下面的这些示例中，A 和 B 数组中长度为 1 的那些轴（缺少的轴会自动补充 1），在广播期间在另一个数组中的相同位置扩展为更大的长度，

A（3d 阵列）：15 x 3 x 5

B（3d 阵列）：15 x 1 x 5

结果（3d 数组）：15 x 3 x 5

A（3d 阵列）：15 x 3 x 5

B（二维数组）：3 x 5

结果（3d 数组）：15 x 3 x 5

A（3d 阵列）：15 x 3 x 5

B（二维数组）：3 x 1

结果（3d 数组）：15 x 3 x 5

A（4d 阵列）：8 x 1 x 6 x 1

B（3d 阵列）：7 x 1 x 5

结果（4d 数组）：8 x 7 x 6 x 5

在下面的例子中，第一个数组的形状为(3,3)c语言实现图像处理，第二个数组的形状为(3,)，等价于(1,3),所以先把第二个数组的形状改成(3,3)，相当于把原来的数组沿0轴复制了两次。

〄 广播机制示意图。 MatA = np.array([[1, 2, 3],[4,5,6],[7,8,9]])

MatB = np.array([1, 2, 3])

MatA + MatB

为了更好地理解这个机制，下面再举几个例子。下图中有三行，分别对应三种广播方式。请检查以下代码。

〄 每一行对应一个广播。 a = np.array([0,10,20,30])

b = np.array([0,1,2])

A = np.stack((a,a,a), 轴=1)

B = np.stack((b,b,b,b))

#对应第一种情况

A + B

#对应第二种情况

A + b

a1 = np.array([[0,10,20,30]]).T

#对应第三种情况

a1 + b

下面的例子不满足广播规则，所以无法进行操作。

A（一维数组）：3

B (1d array): 4 # 最后一个轴长度不兼容

A（二维数组）：4 x 3

B (1d array): 4 # 最后一个轴长度不兼容

A（二维数组）：2 x 1

B (3d array): 8 x 4 x 3 #倒数第二个轴长不兼容

〄 不广播的例子。广播机制总结

2维增减维增

在需要添加轴的位置使用np.newaxis或None。

x = np.arange(6).reshape(2,3)

x, x.shape

x1 = x[:,np.newaxis,:]

x1, x1.形状

# 或

x2 = x[:,None,:]

x2, x2.形状

维数压缩 有时需要去除数组中多余的轴c语言实现图像处理，以减少数组的维数。 numpy.squeeze( )x = np.arange(6).reshape(2,1,3)

y = x.squeeze()

xd = x.__array_interface__[‘data’][0]

yd = y.__array_interface__[‘data’][0]

〄查看内存中数据的地址，验证是否指向同一个内存。 3 数组转置（换轴）x = np.arange(9).reshape(3, 3)

y = np.transpose(x) # or y = x.transpose() or x.T

y = np.transpose(x, [1, 0])

x = np.array([3,2,1,0,4,5,6,7,8,9,10,11,12,13,14,15]).reshape(2, 2, 4)

y1 = np.transpose(x, [1, 0, 2])

请参考下图了解这个 3D 数组在内存中的样子以及它的不同视图。关于这一点，文末附上的进阶章节有详细的解释。

〄 注意轴可以改变，但数据不是固定的。 y2 = np.transpose(x, [2, 0, 1])

#一起编码

x = np.array([3,2,1,0,4,5,6,7,8,9,10,11,12,13,14,15]).reshape(2, 2, 4)

y0 = np.transpose(x, [1, 2, 0])

y1 = np.transpose(x, [1, 0, 2])

y2 = np.transpose(x, [2, 0, 1])

看看改变轴后每个数组的元素是什么样子的。请注意，它们都指向同一条数据。

如何对内存中的相同数据使用不同的轴顺序？其实数据还是那个数据，变化的是每个轴上的步幅。

x.strides, y1.strides, y2.strides

#数据还是一样的

id(x.data), id(y1.data), id(y2.data)

再看一个例子，三维数组有三个轴，注意换轴后每个轴的步长。

x = np.arange(16).reshape(2, 2, 4)

y = x.transpose((1, 0, 2))

两个数组的三个轴对应的步数不同。

换轴后，下标也发生了变化，所以换轴前后相同下标指向的数据是不同的。

〄 轴变了，下标也变了。其实轴的意义主要体现在步长上，所以改变轴在一定意义上就是改变步长。实例

RGB 图像数据

〄 32×32 像素的图像。

看下图，从左到右，对应图像数据在内存中的存储，将一维数组转化为三维数组，替换轴。

那么，为什么要换轴呢？因为不同的包有不同的数据需求，为了使用它们，我们需要根据它们的参数需求来调整数据。

有时候，不需要改变轴，只需改变一个轴上元素的顺序，例如，

# 改变轴上元素的顺序

z = x[…, (3, 2, 1, 0)]

4 通用函数 ufunc 函数 x = np.linspace(0, 2*np.pi, 5)

y, z = np.sin(x), np.cos(x)

# 将结果直接传递给输入x

np.sin(x, x)

性能对比导入时间

导入数学

将 numpy 导入为 np

x = [i for i in range(1000000)]

#math.sin

开始 = time.process_time()

for i, t in enumerate(x):

x[i] = math.sin(t)

math_time = time.process_time()-开始

#numpy.sin

x = np.array(x, dtype=np.float64)

开始 = time.process_time()

np.sin(x, x)

numpy_time = time.process_time()-开始

#比较

math_time、numpy_time、math_time/numpy_time

减少操作

np.add.reduce([1,2,3])

np.add.reduce([[1,2,3],[4,5,6]],axis=1)

np.multiply.reduce([[1,2,3],[4,5,6]],axis=1)

累加操作 np.add.accumulate([1,2,3])

np.add.accumulate([[1,2,3],[4,5,6]],axis=1)

自定义ufunc函数#定义一个python函数

定义 ufunc_diy(x):

c, c0, hc = 0.618, 0.518, 1.0

x = x – int(x)

如果 x >= c:

r = 0.0

elif x

r = x / c0 * hc

其他：

r = (c-x) / (c-c0) * hc

返回 r

x = np.linspace(0, 2, 1000000)

ufunc_diy(x)

开始 = time.process_time()

y1 = np.array([ufunc_diy(t) for t in x])

time_1 = time.process_time()-开始

time_1

np.frompyfunc 函数 ufunc = np.frompyfunc(ufunc_diy, 1, 1)

开始 = time.process_time()

y2 = ufunc(x)

time_2 = time.process_time()-开始

time_2

NumPy 的 C 扩展

本文主要介绍两种扩展方式，

ctypes#ufunc.c

void ufunc_diy(double x, double y, int size) {

双xx,r,c=0.618,c0=0.518,hc=1.0;

for(int i=0;i xx = x[i]-(int)(x[i]);

如果 (xx>=c) r=0.0;

else if (xx else r=(c-xx)/(c-c0)*hc;

y[i]=r;

}

}

”’

#ufunc.py

“”” 包装接受 C 双精度数组的 C 库函数的示例

使用 numpy.ctypeslib 输入。 “””

将 numpy 导入为 np

将 numpy.ctypeslib 导入为 npct

从 ctypes 导入 c_int

array_1d_double = npct.ndpointer(dtype=np.double, ndim=1, flags=’CONTIGUOUS’)

# 加载库，使用 numpy 机制

lib = npct.load_library(“lib_ufunc”, “.”)

#设置返回类型和参数类型

lib.ufunc_diy.restype = 无

lib.ufunc_diy.argtypes = [array_1d_double, array_1d_double, c_int]

def ufunc_diy_func(in_array, out_array):

return lib.ufunc_diy(in_array, out_array, len(in_array))

#编译

# gcc -shared -fPIC -O2 ufunc.c -ldl -o lib_ufunc.so

进口时间

将 numpy 导入为 np

导入 ufunc

开始 = time.process_time()

ufunc.ufunc_diy_func(x, x)

end = time.process_time()

print(“ufunc_diy 时间：”, end-start)

#python test_ufunc.py

# ufunc_diy 时间：0.003 – 0.008

赛通

# ufunc_diy.h

void ufunc_diy(double in_array, double out_array, int size);

#ufunc_diy.c

void ufunc_diy(double x, double y, int size) {

双xx,r,c=0.618,c0=0.518,hc=1.0;

for(int i=0;i xx = x[i]-(int)(x[i]);

如果 (xx>=c) r=0.0;

else if (xx else r=(c-xx)/(c-c0)*hc;

y[i]=r;

}

}

# Cython 支持 NumPy

# 在代码中声明 a = np.array([0,10,20,30])

b = np.array([0,1,2])cimport numpy，使用函数。

#_ufunc_cython.pyx_

“”” 使用将 C 双数组作为输入的 C 函数包装示例

来自 Cython “”” 的 Numpy 声明

# cimport numpy 的 Cython 声明

cimport numpy as np

# 如果你想使用 Cython 的 Numpy-C-API

#（此示例并非绝对必要，但很好的做法）

np.import_array()

# c定义我们的c函数的签名

来自“ufunc_diy.h”的cdef extern：

void ufunc_diy (double in_array, double out_array, int size)

# 创建包装代码，带有 numpy 类型注释

def ufunc_diy_func(np.ndarray[double, ndim=1, mode=”c”] in_array not Noa = np.array([0,10,20,30])

b = np.array([0,1,2])ne,

np.ndarray[double, ndim=1, mode=”c”] out_array not None)：

ufunc_diy(np.PyArray_DATA(in_array),

np.PyArray_DATA(out_array),

in_array.shape[0])

# setup.py

从 distutils.core 导入设置，扩展

导入 numpy

从 Cython.Distutils 导入 build_ext

设置（

cmdclass={‘build_ext’: build_ext},

ext_modules=[Extension(“ufunc_cython”,

sources=[“_ufunc_cython.pyx”, “ufunc_diy.c”],

include_dirs=[numpy.get_include()])],

)

# 或

从 distutils.core 导入设置

导入 numpy

从 Cython.Build 导入 cythonize

设置（

ext_modules=cythonize(“_ufunc_cython.pyx”, annotate=True),

include_dirs=[numpy.get_include()]

)

#编译

python setup.py build_ext –inplace

可以看到多了两个文件，一个是_ufunc_cython.c，一个是ufunc_cython.so（如果是windows就是.pyd）。

c文件是cython将pyx文件解析成ac文件，不依赖平台，而so或pyd文件是c文件编译后的动态链接库，依赖平台。

进口时间

将 numpy 导入为 np

导入 ufunc_cython

开始 = time.process_time()

ufunc_cython.ufunc_diy_func(x, x)

end = time.process_time()

print(“ufunc_diy 时间：”, end-start)

实际运行一下就会知道，当C扩展时，性能有时可以提升数百倍。这样一来，用C语言来实现一些功能是值得的。