NumPy C 代码解释#

狂热主义在于,当你忘记你的目标时,加倍你的努力.— 乔治·桑塔亚那

权威人士是可以告诉您比您真正想知道的更多关于某事的人.— 未知

本页尝试解释一些新代码背后的逻辑.这些解释背后的目的是使某人能够比仅仅盯着代码更容易地理解实现背后的想法.也许通过这种方式,算法可以被更多的人改进,借鉴和/或优化.

内存模型#

ndarray 的一个基本方面是,数组被视为从某个位置开始的内存“块”.这种内存的解释取决于 stride 信息.对于 \(N\) 维数组中的每个维度,一个整数( stride )指示必须跳过多少字节才能到达该维度中的下一个元素.除非您有一个单段数组,否则在遍历数组时必须查阅此 stride 信息.编写接受步长的代码并不难,您只需使用 char* 指针,因为步长以字节为单位.另请记住,步长不必是元素大小的单位倍数.另外,请记住,如果数组的维数是 0(有时称为 rank-0 数组),则 strides 和 dimensions 变量为 NULL .

除了 PyArrayObject 的 strides 和 dimensions 成员中包含的结构信息之外,标志还包含有关如何访问数据的重要信息.特别是,当内存位于根据数据类型数组的合适边界上时,会设置 NPY_ARRAY_ALIGNED 标志.即使您有一个 contiguous 内存块,您也不能仅仅假设可以安全地取消对元素的特定于数据类型的指针的引用.只有设置了 NPY_ARRAY_ALIGNED 标志,这才是安全的操作.在某些平台上它可以工作,但在其他平台上,例如 Solaris,它会导致总线错误.如果您计划写入数组的内存区域,还应确保 NPY_ARRAY_WRITEABLE .也可以获得指向不可写内存区域的指针.有时,当未设置 NPY_ARRAY_WRITEABLE 标志时写入内存区域会很粗鲁.其他时候,它可能导致程序崩溃(例如,作为只读内存映射文件的数据区域).

数据类型封装#

参见

数据类型对象 ( dtype )

datatype 是 ndarray 的一个重要的抽象.操作将查找数据类型,以提供操作数组所需的关键功能.此功能在 PyArray_Descr 结构的 f 成员指向的函数指针列表中提供.通过这种方式,只需在 f 成员中提供带有合适的函数指针的 PyArray_Descr 结构,就可以简单地扩展数据类型的数量.对于内置类型,有一些优化绕过了这种机制,但数据类型抽象的目的是允许添加新的数据类型.

作为一种内置数据类型, void 数据类型允许使用包含一个或多个字段的任意 structured types 作为数组元素. field 只是另一种包含偏移量的数据类型对象,该偏移量指向当前结构化类型.为了支持任意嵌套字段,针对 void 类型实现了数据类型访问的几种递归实现.一种常见的用法是循环遍历字典的元素,并根据存储在给定偏移量处的数据类型对象执行特定操作.这些偏移量可以是任意数字. 因此,必须认识到遇到未对齐数据的可能性,并在必要时加以考虑.

N 维迭代器#

参见

数组迭代

在 NumPy 代码中,一个非常常见的操作是需要迭代通用,跨步,N 维数组的所有元素. 这种通用 N 维循环的操作被抽象成迭代器对象的概念. 要编写 N 维循环,只需从 ndarray 创建一个迭代器对象,使用迭代器对象结构的 dataptr 成员,并在迭代器对象上调用宏 PyArray_ITER_NEXT 以移动到下一个元素. next 元素始终按 C-contiguous 顺序排列. 该宏首先特殊处理 C-contiguous,1-D 和 2-D 情况,这些情况工作非常简单.

对于一般情况,迭代通过跟踪迭代器对象中的坐标计数器列表来工作. 在每次迭代中,最后一个坐标计数器都会增加(从 0 开始). 如果此计数器小于该维度中数组大小减 1(预先计算并存储的值),则计数器会增加,并且 dataptr 成员会按该维度的步幅增加,并且宏结束. 如果到达维度的末尾,则最后一个维度的计数器将重置为零,并且 dataptr 通过减去步幅值乘以该维度中元素数量减 1(也预先计算并存储在迭代器对象的 backstrides 成员中)来移回该维度的开头. 在这种情况下,宏不会结束,但是会减少一个局部维度计数器,以便倒数第二个维度取代最后一个维度发挥的作用,并且先前描述的测试再次在倒数第二个维度上执行. 这样, dataptr 针对任意步幅进行了适当的调整.

coordinates 结构的 PyArrayIterObject 成员维护当前的 N 维计数器,除非底层数组是 C-contiguous,在这种情况下,坐标计数将被绕过. index 的 PyArrayIterObject 成员跟踪迭代器的当前平面索引. 它由 PyArray_ITER_NEXT 宏更新.

广播#

参见

广播

在 NumPy 的祖先 Numeric 中,广播是在 ufuncobject.c 深处埋藏的几行代码中实现的. 在 NumPy 中,广播的概念已被抽象,因此可以在多个位置执行. 广播由函数 PyArray_Broadcast 处理. 此函数需要传入 PyArrayMultiIterObject (或二进制等效项). PyArrayMultiIterObject 跟踪广播的维度数量和每个维度的大小,以及广播结果的总大小. 它还跟踪正在广播的数组的数量以及指向正在广播的每个数组的迭代器的指针.

PyArray_Broadcast 函数接受已经定义的迭代器,并使用它们来确定每个维度中的广播形状(要在广播发生的同时创建迭代器,则使用 PyArray_MultiIterNew 函数).然后,调整迭代器,使每个迭代器都认为它正在迭代具有广播大小的数组.这是通过调整迭代器的维度数量以及每个维度中的 shape 来完成的.这样做是因为迭代器步幅也得到了调整.广播仅调整(或添加)长度为 1 的维度.对于这些维度,步幅变量仅设置为 0,以便该数组的迭代器的数据指针不会随着广播操作在扩展维度上进行而移动.

在 Numeric 中,广播始终使用扩展维度的 0 值步幅来实现.在 NumPy 中,它的实现方式完全相同.最大的区别在于,步幅数组现在保存在 PyArrayIterObject 中,广播结果中涉及的迭代器保存在 PyArrayMultiIterObject 中,并且 PyArray_Broadcast 调用实现了 通用广播规则 .

数组标量#

参见

标量

数组标量提供了一个 Python 类型层次结构,允许数组中存储的数据类型与从数组中提取元素时返回的 Python 类型之间存在一一对应的关系.对象数组是此规则的一个例外.对象数组是任意 Python 对象的异构集合.当您从对象数组中选择一个项目时,您会得到原始的 Python 对象(而不是对象数组标量,它确实存在,但在实际应用中很少使用).

数组标量还提供与数组相同的方法和属性,目的是可以使用相同的代码来支持任意维度(包括 0-维度).数组标量是只读的(不可变的),但 void 标量除外,它也可以写入,以便结构化数组字段设置可以更自然地工作( a[0]['f1'] = value ).

索引#

参见

在 ndarrays 上的索引 , 索引例程

所有 Python 索引操作 arr[index] 首先通过准备索引并找到索引类型来组织.支持的索引类型有:

  • 整数

  • newaxis

  • slice

  • Ellipsis

  • 整数数组/类数组(高级)

  • 布尔值(单个布尔数组);如果作为索引的布尔数组多于一个,或者形状不完全匹配,则布尔数组将转换为整数数组.

  • 0-d 布尔值(以及整数);0-d 布尔数组是一个特殊情况,必须在高级索引代码中处理.它们表示 0-d 布尔数组必须解释为整数数组.

以及标量数组的特殊情况,表示整数数组被解释为整数索引,这很重要,因为整数数组索引强制复制,但如果返回标量则会被忽略(完整整数索引).准备好的索引保证有效,但超出范围的值和高级索引的广播错误除外.这包括为不完整的索引添加 Ellipsis ,例如当使用单个整数索引二维数组时.

下一步取决于找到的索引的类型.如果所有维度都用整数索引,则返回或设置标量.单个布尔索引数组将调用专门的布尔函数.包含 Ellipsis 或 slice 但没有高级索引的索引将始终通过计算新的步幅和内存偏移来创建旧数组的视图.然后,可以使用 PyArray_CopyObject 返回此视图或对于赋值进行填充.请注意,在其他分支中也可能在临时数组上调用 PyArray_CopyObject ,以支持当数组是对象 dtype 时复杂的赋值.

高级索引#

到目前为止,最复杂的情况是高级索引,它可能与基于视图的典型索引结合使用,也可能不结合使用. 这里,整数索引被解释为基于视图的索引. 在尝试理解这一点之前,您可能需要熟悉它的微妙之处. 高级索引代码有三个不同的分支和一个特殊情况:

  • 有一个索引数组,它可以和赋值数组一样被简单地迭代. 例如,它们可能是连续的. 此外,索引数组必须是 intp 类型,赋值中的值数组应该是正确的类型. 这纯粹是一个快速通道.

  • 只有整数数组索引,因此不存在子数组.

  • 混合了基于视图的索引和高级索引. 在这种情况下,基于视图的索引定义了一组子数组,这些子数组通过高级索引组合在一起. 例如, arr[[1, 2, 3], :] 是通过垂直堆叠子数组 arr[1, :] , arr[2, :] 和 arr[3, :] 创建的.

  • 存在一个子数组,但它只有一个元素. 这种情况可以被视为没有子数组,但在设置期间需要小心.

决定适用哪种情况,检查广播以及确定所需的转置类型都在 PyArray_MapIterNew 中完成. 设置完成后,有两种情况. 如果没有子数组或只有一个元素,则不需要子数组迭代,并且准备一个迭代器,该迭代器迭代所有索引数组以及结果或值数组. 如果存在子数组,则准备三个迭代器. 一个用于索引数组,一个用于结果或值数组(减去其子数组),另一个用于原始数组和结果/赋值数组的子数组. 前两个迭代器给出(或允许计算)指向子数组开始位置的指针,然后可以重新启动子数组迭代.

当高级索引彼此相邻时,可能需要转置. 所有必要的转置都由 PyArray_MapIterSwapAxes 处理,并且必须由调用者处理,除非要求 PyArray_MapIterNew 分配结果.

准备之后,获取和设置相对简单,尽管需要考虑不同的迭代模式. 除非在项目获取期间只有一个索引数组,否则会预先检查索引的有效性. 否则,为了优化,它会在内部循环本身中处理.

通用函数#

参见

通用函数 ( ufunc ) , 通用函数( ufunc )基础

通用函数是可调用对象,它接受 \(N\) 个输入并产生 \(M\) 个输出,方法是将基本的 1-D 循环(逐个元素工作)包装成完整的,易于使用的函数,这些函数无缝地实现 broadcasting , type-checking , buffered coercion 和 output-argument handling . 新的通用函数通常在 C 中创建,尽管有一种机制可以从 Python 函数创建 ufunc( frompyfunc ). 用户必须提供一个 1-D 循环,该循环实现基本函数,该函数接受输入标量值并将结果标量放入适当的输出槽中,如实现中所述.

设置#

每个 ufunc 计算都涉及与设置计算相关的一些开销. 这种开销的实际意义是,即使 ufunc 的实际计算非常快,您也可以编写数组和类型特定的代码,这些代码对于小数组会比 ufunc 更快. 特别是,使用 ufunc 对 0-D 数组执行许多计算将比其他基于 Python 的解决方案慢(静默导入的 scalarmath 模块的存在正是为了使数组标量具有基于 ufunc 的计算的外观和感觉,同时显著减少开销).

当调用一个 ufunc 时,有很多事情必须完成.从这些设置操作中收集的信息存储在一个循环对象中.该循环对象是一个 C 结构(可以成为一个 Python 对象,但没有初始化为 Python 对象,因为它只在内部使用).这个循环对象具有与 PyArray_Broadcast 一起使用的布局,因此广播可以以与其他代码段相同的方式处理.

首先要做的是在线程特定的全局字典中查找缓冲区大小,错误掩码和相关错误对象的当前值.错误掩码的状态控制发现错误情况时会发生什么.应该注意的是,硬件错误标志的检查仅在每次执行一维循环后才进行.这意味着如果输入和输出数组是连续的并且类型正确,以便执行单个一维循环,则可能要等到计算完数组的所有元素后才检查这些标志.在线程特定的字典中查找这些值会花费时间,对于所有但非常小的数组来说,这些时间很容易被忽略.

在检查线程特定的全局变量之后,对输入进行评估以确定 ufunc 应该如何进行,并在必要时构造输入和输出数组.任何不是数组的输入都将转换为数组(如果需要,使用上下文).记录哪些输入是标量(因此转换为 0-D 数组).

接下来,根据输入数组类型,从 ufunc 可用的 1-D 循环中选择合适的 1-D 循环.选择此 1-D 循环的方法是尝试将输入的数据类型签名与可用签名进行匹配.与内置类型对应的签名存储在 ufunc 结构的 ufunc.types 成员中.与用户定义的类型对应的签名存储在函数信息的链表中,头元素存储为 CObject ,位于以数据类型编号为键的 userloops 字典中(参数列表中的第一个用户定义类型用作键).搜索签名,直到找到一个签名,输入数组都可以安全地强制转换为该签名(忽略任何不允许确定结果类型的标量参数).这种搜索过程的含义是,在存储签名时,“较小的类型”应放置在“较大的类型”之下.如果未找到 1-D 循环,则报告错误.否则, argument_list 将使用存储的签名进行更新—以防需要强制转换并修复 1-D 循环假定的输出类型.

如果 ufunc 有 2 个输入和 1 个输出,并且第二个输入是 Object 数组,则会执行一个特殊情况检查,以便如果第二个输入不是 ndarray,具有 __array_priority__ 属性,并且具有 __r{op}__ 特殊方法,则返回 NotImplemented .这样,Python 会收到信号,让另一个对象有机会完成操作,而不是使用通用对象数组计算.这允许(例如)稀疏矩阵覆盖乘法运算符 1-D 循环.

对于小于指定缓冲区大小的输入数组,将复制所有非连续,未对齐或字节序错误的数组,以确保对于小数组,使用单个循环.然后,为所有输入数组创建数组迭代器,并将生成的迭代器集合广播到单个形状.

然后处理输出参数(如果有),并构造任何缺失的返回数组.如果任何提供的输出数组没有正确的类型(或未对齐)并且小于缓冲区大小,则将构造一个新的输出数组,并设置特殊的 NPY_ARRAY_WRITEBACKIFCOPY 标志.在该函数的末尾,调用 PyArray_ResolveWritebackIfCopy 以便将其内容复制回输出数组.然后处理输出参数的迭代器.

最后,确定了如何执行循环机制,以确保所有输入数组的元素组合在一起,以生成正确类型的输出数组. 循环执行的选项包括单循环(适用于 contiguous ,对齐且数据类型正确的情况),跨步循环(适用于非连续但仍对齐且数据类型正确的情况)和缓冲循环(适用于未对齐或数据类型错误的情况). 根据所需的执行方法,然后设置和计算循环.

函数调用#

本节介绍如何为三种不同的执行方式设置和执行基本的通用函数计算循环. 如果在编译期间定义了 NPY_ALLOW_THREADS ,那么只要不涉及对象数组,则在调用循环之前,Python 全局解释器锁 (GIL) 将被释放. 如果需要处理错误情况,则会重新获取它. 硬件错误标志仅在 1-D 循环完成后才进行检查.

单循环#

这是所有情况中最简单的一种. 通过完全调用一次底层的 1-D 循环来执行 ufunc. 只有当输入和输出的对齐数据类型正确(包括字节顺序)且所有数组都具有统一步幅( contiguous ,0-D 或 1-D)时,才有可能. 在这种情况下,将调用一次 1-D 计算循环,以计算整个数组的计算. 请注意,仅在整个计算完成后才检查硬件错误标志.

跨步循环#

当输入和输出数组对齐且类型正确,但步幅不统一(非连续且 2-D 或更大)时,则采用第二个循环结构进行计算. 这种方法会将输入和输出参数的所有迭代器转换为迭代除最大维度之外的所有维度. 然后,内部循环由底层的 1-D 计算循环处理. 外部循环是转换后的迭代器的标准迭代器循环. 在每个 1-D 循环完成后检查硬件错误标志.

缓冲循环#

每当输入和/或输出数组与底层 1-D 循环期望的值未对齐或数据类型错误(包括字节交换)时,此代码将处理这种情况. 数组也被假定为非连续的. 该代码的工作方式与跨步循环非常相似,不同之处在于内部 1-D 循环被修改,以便对输入执行预处理,并对输出执行 bufsize 块的后处理(其中 bufsize 是用户可设置的参数). 在复制过来的数据上调用底层的 1-D 计算循环(如果需要). 在这种情况下,设置代码和循环代码要复杂得多,因为它必须处理:

  • 临时缓冲区的内存分配

  • 决定是否在输入和输出数据上使用缓冲区(未对齐和/或数据类型错误)

  • 复制和可能转换任何需要缓冲区的输入或输出的数据.

  • 特殊情况 Object 数组,以便在需要复制和/或强制转换时正确处理引用计数.

  • 将内部 1-D 循环分解为 bufsize 块(可能存在余数).

同样,在每个 1-D 循环结束时检查硬件错误标志.

最终输出操作#

Ufuncs 允许其他类似数组的类无缝地通过接口传递,因为特定类的输入将导致输出与该类相同. 其工作原理如下. 如果任何输入不是 ndarrays 并且定义了 __array_wrap__ 方法,则具有最大 __array_priority__ 属性的类将确定所有输出的类型(除了传入的任何输出数组). 输入数组的 __array_wrap__ 方法将被调用,并将从 ufunc 返回的 ndarray 作为其输入. 支持两种调用 __array_wrap__ 函数的样式. 第一种将 ndarray 作为第一个参数,将“上下文”的元组作为第二个参数. 上下文是 (ufunc, arguments, output argument number). 这是第一次尝试调用. 如果发生 TypeError ,则仅使用 ndarray 作为第一个参数来调用该函数.

方法#

有三种 ufunc 方法需要进行与通用 ufunc 类似的计算.它们是 ufunc.reduce , ufunc.accumulate 和 ufunc.reduceat .这些方法中的每一种都需要一个设置命令,然后是一个循环.这些方法有四种可能的循环样式,分别对应于无元素,单元素,跨步循环和缓冲循环.这些是与通用函数调用中实现的相同的基本循环样式,除了无元素和单元素的情况,这两种情况是分别在输入数组对象具有 0 和 1 个元素时发生的特殊情况.

设置#

所有三种方法的设置函数都是 construct_reduce .此函数创建一个归约循环对象,并用完成循环所需的参数填充它.所有这些方法仅适用于采用 2 个输入并返回 1 个输出的 ufunc.因此,底层 1-D 循环被选择,假设签名为 [otype, otype, otype] ,其中 otype 是请求的归约数据类型.然后,从(每个线程)全局存储中检索缓冲区大小和错误处理.对于未对齐或具有不正确数据类型的小数组,会进行复制,以便使用代码的非缓冲部分.然后,选择循环策略.如果数组中有 1 个元素或 0 个元素,则选择简单的循环方法.如果数组未对齐且具有正确的数据类型,则选择跨步循环.否则,必须执行缓冲循环.然后建立循环参数,并构造返回数组.输出数组的 shape 不同,具体取决于方法是 reduce , accumulate 还是 reduceat .如果已经提供了输出数组,则检查其形状.如果输出数组不是 C-contiguous,对齐的且具有正确的数据类型,则会使用 NPY_ARRAY_WRITEBACKIFCOPY 标志创建一个临时副本.通过这种方式,这些方法将能够使用行为良好的输出数组,但当在函数完成时调用 PyArray_ResolveWritebackIfCopy 时,结果将被复制回真正的输出数组中.最后,设置迭代器以循环遍历正确的 axis (取决于提供给该方法的轴的值),并且设置例程返回到实际的计算例程.

Reduce#

所有 ufunc 方法都使用相同的底层 1-D 计算循环,输入和输出参数经过调整,以便进行适当的归约. 例如, reduce 起作用的关键是,调用 1-D 循环时,输出和第二个输入指向内存中的同一位置,并且两者的步长均为 0.第一个输入指向输入数组,其步长由为所选轴提供的适当跨度给出. 这样,执行的操作是

\begin{align} o & = & i[0] \\ o & = & i[k]\textrm{<op>}o\quad k=1\ldots N \end{align}

其中 \(N+1\) 是输入 \(i\) 中的元素数, \(o\) 是输出, \(i[k]\) 是沿所选轴的 \(k^{\textrm{th}}\) 的第 \(i\) 个元素. 对于具有大于 1 维的数组,此基本操作会重复执行,以便沿所选轴为每个 1-D 子数组进行归约. 删除所选维度的迭代器处理此循环.

对于缓冲循环,必须注意在调用循环函数之前复制和转换数据,因为底层循环需要已对齐且具有正确数据类型(包括字节顺序)的数据. 缓冲循环必须在调用循环函数之前,处理不大于用户指定的 bufsize 的块上的复制和强制转换.

Accumulate#

accumulate 方法与 reduce 方法非常相似,因为输出和第二个输入都指向输出.不同之处在于,第二个输入指向当前输出指针之后一个步长的内存.因此,执行的操作是

\begin{align} o[0] & = & i[0] \\ o[k] & = & i[k]\textrm{<op>}o[k-1]\quad k=1\ldots N. \end{align}

输出的形状与输入相同,当所选轴的形状为 \(N\) 时,每个 1-D 循环对 \(N+1\) 个元素进行运算.同样,缓冲循环会在调用底层 1-D 计算循环之前,注意复制和转换数据.

Reduceat#

reduceat 函数是 reduce 和 accumulate 函数的概括.它实现了由索引指定的输入数组范围上的 reduce .在进行循环计算之前,会检查额外的 indices 参数,以确保每个输入对于所选维度上的输入数组而言不是太大.循环实现的处理使用的代码与 reduce 代码非常相似,该代码重复的次数与 indices 输入中的元素数一样多.特别是:传递给底层 1-D 计算循环的第一个输入指针指向由索引数组指示的输入数组中的正确位置.此外,传递给底层 1-D 循环的输出指针和第二个输入指针指向内存中的相同位置.1-D 计算循环的大小固定为当前索引和下一个索引之间的差值(当当前索引是最后一个索引时,则下一个索引被假定为数组沿所选维度的长度).通过这种方式,1-D 循环将实现指定索引上的 reduce .

未对齐或循环数据类型与输入和/或输出数据类型不匹配的情况,会使用缓冲代码进行处理,其中数据被复制到临时缓冲区,并在调用底层 1-D 函数之前根据需要转换为正确的数据类型.临时缓冲区以(元素)大小创建,不大于用户可设置的 buffer-size 值.因此,循环必须足够灵活,能够多次调用底层 1-D 计算循环,以便在不大于 buffer-size 的块中完成总计算.