NumPy C 代码解释#

狂热就是当你遗忘了目标时，反而加倍努力。 — 乔治·桑塔亚纳 (George Santayana)

所谓权威，就是能告诉你关于某件事比你真正想知道的还要多的人。 — 佚名 (Unknown)

本页旨在解释一些新代码背后的逻辑。这些解释的目的在于使人们能够比仅盯着代码更容易地理解实现背后的思路。也许通过这种方式，这些算法可以被更多的人改进、借鉴和/或优化。

内存模型#

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

除了 PyArrayObject 的步长 (strides) 和维度 (dimensions) 成员中包含的结构信息之外，标志 (flags) 还包含有关如何访问数据的重要信息。特别地，当内存处于适合数据类型数组的边界时，会设置 NPY_ARRAY_ALIGNED 标志。即使您有一个 连续的 (contiguous) 内存块，也不能简单地假定解引用指向该元素的特定数据类型指针是安全的。只有当设置了 NPY_ARRAY_ALIGNED 标志时，这才是安全的操作。在某些平台上它可以工作，但在其他平台（如 Solaris）上，它将导致总线错误。如果您计划写入数组的内存区域，还应确保设置了 NPY_ARRAY_WRITEABLE 标志。获取指向不可写内存区域的指针也是可能的。有时，在未设置 NPY_ARRAY_WRITEABLE 标志时写入该内存区域只会被视作粗鲁的操作。而有时它会导致程序崩溃（例如：作为只读内存映射文件的核心数据区）。

数据类型封装#

另请参阅

数据类型对象 (dtype)

数据类型 (datatype) 是 ndarray 的重要抽象。操作将依靠数据类型来提供操作数组所需的核心功能。该功能由 PyArray_Descr 结构的 f 成员所指向的函数指针列表提供。通过这种方式，只需提供一个在 f 成员中具有合适函数指针的 PyArray_Descr 结构，就可以扩展数据类型的数量。对于内置类型，有一些优化会绕过此机制，但数据类型抽象的重点是允许添加新的数据类型。

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

N 维迭代器#

另请参阅

遍历数组

在大部分 NumPy 代码中，一个非常常见的操作是需要迭代一个通用的、跨步的 (strided) N 维数组的所有元素。这种通用 N 维循环的操作在迭代器对象概念中得到了抽象。要编写 N 维循环，您只需从 ndarray 创建一个迭代器对象，操作该迭代器对象结构的 dataptr 成员，并在迭代器对象上调用宏 PyArray_ITER_NEXT 以移动到下一个元素。下一个 (next) 元素始终采用 C 连续顺序。该宏的工作原理是：首先将 C 连续、1 维和 2 维情况作为特例处理，这些情况的处理非常简单。

对于一般情况，迭代通过在迭代器对象中跟踪一个坐标计数器列表来工作。在每次迭代中，最后一个坐标计数器增加（从 0 开始）。如果该计数器小于该维度数组大小减一（一个预先计算并存储的值），则增加该计数器，并且 dataptr 成员增加该维度上的步长，宏结束。如果到达维度的末尾，则将最后一个维度的计数器重置为零，并且通过减去步长值乘以该维度元素数减一，将 dataptr 移回该维度的起始位置（这同样是预先计算好的，并存储在迭代器对象的 backstrides 成员中）。在这种情况下，宏并不结束，而是递减本地维度计数器，以便倒数第二个维度代替最后一个维度的角色，并在倒数第二个维度上再次执行前面描述的测试。通过这种方式，可以针对任意跨步对 dataptr 进行适当调整。

PyArrayIterObject 结构的 coordinates 成员维护当前的 N 维计数器，除非底层数组是 C 连续的，在这种情况下会绕过坐标计数。PyArrayIterObject 的 index 成员跟踪迭代器的当前扁平索引。它由 PyArray_ITER_NEXT 宏更新。

广播#

另请参阅

广播

在 NumPy 的前身 Numeric 中，广播 (broadcasting) 仅是用深埋在 ufuncobject.c 中的几行代码实现的。在 NumPy 中，广播的概念已经抽象出来，以便可以在多个地方执行。广播由函数 PyArray_Broadcast 处理。该函数需要传入一个 PyArrayMultiIterObject（或二进制等效结构）。PyArrayMultiIterObject 跟踪广播后的维度数和每个维度的大小，以及广播结果的总大小。它还跟踪正在进行广播的数组数量，以及指向每个广播数组的迭代器的指针。

PyArray_Broadcast 函数获取已定义的迭代器，并使用它们来确定每个维度上的广播形状（要在广播发生的同时创建迭代器，则使用 PyArray_MultiIterNew 函数）。然后，对迭代器进行调整，使每个迭代器都认为自己正在迭代具有广播大小的数组。这是通过调整迭代器的维度数量以及每个维度的 形状 (shape) 来实现的。这之所以可行，是因为迭代器步长也进行了调整。广播仅调整（或添加）长度为 1 的维度。对于这些维度，步长变量仅设置为 0，以便该数组的迭代器的数据指针在广播操作在扩展维度上运行时不会移动。

在 Numeric 中，广播总是通过对扩展维度使用 0 值步长来实现的。在 NumPy 中也是以完全相同的方式完成的。最大的区别在于，现在步长数组在 PyArrayIterObject 中被跟踪，参与广播结果的迭代器在 PyArrayMultiIterObject 中被跟踪，并且 PyArray_Broadcast 调用实现了 一般广播规则。

数组标量#

另请参阅

标量

数组标量 (array scalars) 提供了一个 Python 类型层次结构，该结构允许在数组中存储的数据类型与从数组中提取元素时返回的 Python 类型之间建立一一对应的关系。对象数组 (object arrays) 是该规则的一个例外。对象数组是任意 Python 对象的异构集合。当您从对象数组中选择一个项时，您将获得原始 Python 对象（而不是对象数组标量，尽管它确实存在，但很少用于实际目的）。

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

索引#

另请参阅

ndarray 上的索引，索引例程

所有 Python 索引操作 arr[index] 都是通过首先准备索引并确定索引类型来组织的。支持的索引类型包括

  • 整数 (integer)

  • newaxis (新轴)

  • 切片 (slice)

  • 省略号 (Ellipsis)

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

  • 布尔（单个布尔数组）；如果索引中包含多个布尔数组，或者形状不完全匹配，则布尔数组将被转换为整数数组。

  • 0 维布尔（以及整数）；0 维布尔数组是一个特殊情况，必须在高级索引代码中处理。它们表示 0 维布尔数组必须被解释为整数数组。

以及标量数组特例，表示整数数组被解释为整数索引，这很重要，因为整数数组索引会强制进行复制，但如果返回的是标量（全整数索引），则会被忽略。除了高级索引的越界值和广播错误之外，准备好的索引保证是有效的。例如，当用单个整数索引二维数组时，这包括为不完整的索引添加一个 Ellipsis。

下一步取决于所发现的索引类型。如果所有维度都用整数索引，则返回或设置一个标量。单个布尔索引数组将调用专门的布尔函数。包含 Ellipsis 或 切片 (slice) 但不包含高级索引的索引，将始终通过计算新的步长和内存偏移量来创建旧数组的视图。然后，该视图要么被返回，要么（对于赋值操作）使用 PyArray_CopyObject 进行填充。注意，在其他分支中也可以在临时数组上调用 PyArray_CopyObject，以便在数组为对象 dtype 时支持复杂的赋值操作。

高级索引#

到目前为止，最复杂的情况是高级索引，它可能会也可能不会与典型的基于视图的索引相结合。在这里，整数索引被解释为基于视图的。在尝试理解这一点之前，您可能希望使自己熟悉其细节。高级索引代码有三个不同的分支和一个特殊情况

  • 存在一个索引数组，它以及赋值数组都可以进行简单迭代。例如，它们可以是连续的。此外，索引数组必须是 intp 类型，并且赋值中的值数组应当是正确的类型。这完全是一个快速路径 (fast path)。

  • 仅存在整数数组索引，因此不存在子数组。

  • 基于视图的索引与高级索引混合。在这种情况下，基于视图的索引定义了由高级索引合并而成的子数组集合。例如，arr[[1, 2, 3], :] 是通过垂直堆叠子数组 arr[1, :]、arr[2, :] 和 arr[3, :] 来创建的。

  • 存在一个子数组，但它恰好只有一个元素。这种情况可以按照不存在子数组的方式进行处理，但在设置期间需要注意。

判定适用哪种情况、检查广播以及确定所需的转置类型，都在 PyArray_MapIterNew 中完成。设置完成后，存在两种情况。如果没有子数组或者子数组仅有一个元素，则不需要子数组迭代，并准备一个迭代器，该迭代器迭代所有索引数组 以及 结果或值数组。如果有子数组，则准备三个迭代器。一个用于索引数组，一个用于结果或值数组（减去其子数组），还有一个用于原始数组和结果/赋值数组的子数组。前两个迭代器给出了（或允许计算出）指向子数组起始位置的指针，然后就可以重新开始子数组迭代。

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

准备就绪后，获取和设置相对直接，尽管需要考虑不同的迭代模式。除非在获取项期间只有一个索引数组，否则会预先检查索引的有效性。否则，为了优化，它将在内层循环本身中进行处理。

通用函数#

另请参阅

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

通用函数是可调用对象，它们接受 \(N\) 个输入并产生 \(M\) 个输出。它们通过将逐个元素进行计算的基本 1 维循环包装成易于使用的完整函数，从而无缝地实现 广播、类型检查、缓冲强制类型转换 以及 输出参数处理。新的通用函数通常是用 C 语言创建的，尽管也存在一种从 Python 函数创建 ufunc 的机制（frompyfunc）。正如在实现中所说明的那样，用户必须提供一个实现基本函数的 1 维循环，该循环接收输入标量值并将产生的结果标量放入相应的输出槽中。

设置 (Setup)#

每个 ufunc 计算都涉及一些与设置计算相关的开销。这种开销的实际意义在于，即使 ufunc 的实际计算非常快，您也能编写针对特定数组和类型的代码，对于小数组来说，这些代码运行起来比 ufunc 更快。特别是，使用 ufunc 对 0 维数组执行许多计算将比其他基于 Python 的解决方案更慢（静默导入的 scalarmath 模块的存在恰恰是为了给数组标量提供基于 ufunc 的计算的外观和感觉，同时显著减少开销）。

当调用 ufunc 时，必须做许多事情。从这些设置操作中收集的信息存储在一个循环对象中。这个循环对象是一个 C 结构（它可以成为 Python 对象，但由于它仅在内部使用，因此并未以此进行初始化）。该循环对象具有与 PyArray_Broadcast 一起使用所需的布局，以便能以与其他代码部分相同的方式处理广播。

首先在特定于线程的全局字典中查找缓冲区大小、错误掩码以及相关错误对象的当前值。错误掩码的状态控制发生错误情况时的处理。需要注意的是，对硬件错误标志的检查仅在执行完每个 1 维循环后进行。这意味着如果输入和输出数组是连续的且类型正确，从而执行了单个 1 维循环，那么在计算完数组的所有元素之前可能不会检查标志。在线程特定的字典中查找这些值需要时间，对于除极小数组外的所有数组，这一点都很容易被忽略。

检查线程特定的全局变量后，将对输入进行评估以确定 ufunc 应如何继续，并在必要时构建输入和输出数组。任何不是数组的输入都将被转换为数组（必要时使用上下文）。记录哪些输入是标量（并因此被转换为 0 维数组）。

接着，根据输入数组类型，从 ufunc 可用的 1 维循环中选择一个合适的 1 维循环。通过尝试将输入数据类型的签名与可用签名进行匹配来选择此 1 维循环。与内置类型对应的签名存储在 ufunc 结构的 ufunc.types 成员中。与用户自定义类型对应的签名存储在函数信息的链表中，其头部元素作为 CObject 存储在 userloops 字典中，其键为数据类型编号（参数列表中第一个用户自定义的类型用作键）。搜索这些签名，直到找到一个输入数组都可以安全转型的签名（忽略不允许决定结果类型的任何标量参数）。此搜索过程的含义是，在存储签名时，“较小类型”应放在“较大类型”下方。如果未找到 1 维循环，则报告错误。否则，使用存储的签名更新 argument_list — 以备转型需要并确定 1 维循环所假定的输出类型。

如果 ufunc 具有 2 个输入和 1 个输出，且第二个输入是 Object 数组，则会执行特殊情况检查：如果第二个输入不是 ndarray、具有 __array_priority__ 属性且具有 __r{op}__ 特殊方法，则返回 NotImplemented。通过这种方式通知 Python，给另一个对象一个完成操作的机会，而不是使用通用的对象数组计算。这允许（例如）稀疏矩阵覆盖乘法运算符的 1 维循环。

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

然后处理输出参数（如果有），并构建任何缺失的返回数组。如果任何提供的输出数组类型不正确（或未对齐）且小于缓冲区大小，则构建一个新的输出数组，并设置特殊的 NPY_ARRAY_WRITEBACKIFCOPY 标志。在函数结束时，将调用 PyArray_ResolveWritebackIfCopy，以便将其内容复制回输出数组中。接着处理输出参数的迭代器。

最后，决定如何执行循环机制，以确保将输入数组的所有元素合并来产生正确类型的输出数组。循环执行的选项有：单循环（用于连续、对齐且数据类型正确的情况）、跨步循环（用于非连续但仍然对齐且数据类型正确的情况）以及缓冲循环（用于未对齐或数据类型不正确的情况）。根据所需的执行方法，随后将设置并计算循环。

函数调用#

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

单循环#

这是所有情况中最简单的一种。通过精确调用一次底层 1 维循环来执行 ufunc。只有在输入和输出的数据都是对齐的且类型正确（包括字节序），并且所有数组都具有均匀的步长（连续、0 维或 1 维）时，这才可能。在这种情况下，1 维计算循环会被调用一次，以计算整个数组。注意，只有在整个计算完成后才会检查硬件错误标志。

跨步循环#

当输入和输出数组对齐且类型正确，但跨步不均匀（非连续且为 2 维或更大维）时，计算将采用第二种循环结构。该方法转换输入和输出参数的所有迭代器，以便对除最大维度之外的所有维度进行迭代。内层循环随后由底层的 1 维计算循环处理。外层循环是针对转换后的迭代器的标准迭代器循环。在每个 1 维循环完成后，会检查硬件错误标志。

缓冲循环#

这段代码用于处理输入和/或输出数组未对齐，或者与底层 1 维循环期望的数据类型不匹配（包括字节交换过）的情况。这些数组也被假定为非连续的。其代码运行方式与跨步循环非常相似，不同之处在于对内层 1 维循环进行了修改，从而以大小为 bufsize 的块（其中 bufsize 是用户可设置的参数）对输入进行预处理，并对输出进行后处理。在被复制过来的数据上（如有需要）调用底层的 1 维计算循环。在这种情况下，设置代码和循环代码要复杂得多，因为它必须处理：

  • 临时缓冲区的内存分配

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

  • 为任何需要缓冲区的输入或输出复制并可能转型数据。

  • 对 Object 数组进行特殊处理，以便在需要进行复制和/或转型时正确处理引用计数。

  • 将内层 1 维循环分解为 bufsize 大小的块（可能有余数）。

同样，在每个 1 维循环结束时都会检查硬件错误标志。

最终输出操作#

Ufunc 允许其他类似数组的类无缝通过该接口，因为特定类的输入将导致输出也属于该类。其工作机制如下。如果任何输入不是 ndarray 并且定义了 __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。这些方法中的每一个都需要一条设置命令，后跟一个循环。这些方法可能对应四种循环样式：无元素、单元素、跨步循环和缓冲循环。除了无元素 and 单元素情况外，这些循环样式与针对通用函数调用所实现的循环样式基本相同；无元素和单元素情况是分别在输入数组对象具有 0 个 and 1 个元素时发生的特殊情况。

设置 (Setup)#

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

Reduce#

所有的 ufunc 方法都使用相同的底层 1 维计算循环，并调整了输入和输出参数，以便进行适当的规约。例如，reduce 运行的关键在于：调用 1 维循环时，输出和第二个输入指向内存中的相同位置，且两者的步长都为 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]\) 是 \(i\) 沿所选轴的第 \(k\) 个元素。对于大于 1 维的数组，会重复此基本操作，以便对沿所选轴的每个 1 维子数组进行规约。移除了所选维度的迭代器将处理此循环。

对于缓冲循环，在调用循环函数之前必须注意复制和转型数据，因为底层循环需要对齐且数据类型正确（包括字节序）的数据。在对不大于用户指定 bufsize 的块调用循环函数之前，缓冲循环必须处理好这种复制和转型。

Accumulate#

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

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

输出与输入具有相同的形状，并且当所选轴上的形状为 \(N+1\) 时，每个 1 维循环都会对 \(N\) 个元素进行操作。同样，在调用底层 1 维计算循环之前，缓冲循环会注意复制并转型数据。

Reduceat#

reduceat 函数是 reduce 和 accumulate 函数的推广。它在索引指定的输入数组范围上实现 reduce。在进行循环计算之前，会对额外的 indices 参数进行检查，以确保每个输入对于沿所选维度的输入数组来说都不会太大。循环实现使用与 reduce 非常相似的代码进行处理，并重复 indices 输入中的元素个数相同的次数。特别地：传给底层 1 维计算循环的第一个输入指针指向索引数组指示的输入数组的正确位置。此外，传给底层 1 维循环的输出指针和第二个输入指针指向内存中的相同位置。1 维计算循环的大小固定为当前索引与下一个索引之间的差值（当前索引是最后一个索引时，下一个索引被假定为数组沿所选维度的长度）。通过这种方式，1 维循环将在指定的索引上实现 reduce。

未对齐或循环数据类型与输入和/或输出数据类型不匹配的情况，将使用缓冲代码来处理，其中在调用底层 1 维函数之前，数据会被复制到临时缓冲区中，并在必要时转换为正确的数据类型。创建的临时缓冲区大小（元素个数）不大于用户可设置的缓冲区大小值。因此，该循环必须具有足够的灵活性，以便在不大于缓冲区大小的块中调用足够次数的底层 1 维计算循环来完成总的计算。