NumPy C 代码解释#

狂热在于当你忘记了自己的目标时加倍努力。 — 乔治·桑塔亚纳

权威是指一个比你更了解某件事的人，而你并不想了解那么详细。 — 未知

本页试图解释一些新代码段背后的逻辑。这些解释的目的是让人们能够更容易地理解实现背后的想法，而不是仅仅盯着代码看。也许通过这种方式，算法可以被更多的人改进、借用和/或优化。

内存模型#

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

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

数据类型封装#

另请参阅

数据类型对象 (dtype)

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

内置数据类型之一，void 数据类型允许将 1 个或多个字段作为数组的元素的任意结构化类型。一个字段 仅仅是另一个数据类型对象以及到当前结构化类型中的偏移量。为了支持任意嵌套字段，为 void 类型实现了若干个数据类型访问的递归实现。一个常见的习惯用法是遍历字典中的元素，并根据存储在给定偏移量处的对象类型对象执行特定操作。这些偏移量可以是任意数字。因此，如果需要，必须识别并考虑遇到未对齐数据的可能性。

N 维迭代器#

另请参阅

遍历数组

在许多 NumPy 代码中，一个非常常见的操作是需要遍历一个通用、有步长、N 维数组的所有元素。通用 N 维循环的操作在迭代器对象的概念中得到了抽象。要编写一个 N 维循环，你只需从 ndarray 创建一个迭代器对象，使用迭代器对象结构的dataptr 成员，并调用宏PyArray_ITER_NEXT 来移动到下一个元素。next 元素始终以 C 连续顺序排列。该宏的工作原理是首先对 C 连续、1 维和 2 维情况进行特殊处理，这些情况非常简单。

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

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

广播#

另请参阅

广播

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

PyArray_Broadcast 函数使用已定义的迭代器来确定每个维度中的广播形状（要同时创建迭代器和广播，然后使用PyArray_MultiIterNew 函数）。然后，迭代器会进行调整，以便每个迭代器认为它正在遍历一个具有广播大小的数组。这是通过调整迭代器的维度数量和每个维度中的形状 来完成的。这之所以有效是因为迭代器的步长也进行了调整。广播只调整（或添加）长度为 1 的维度。对于这些维度，步长变量 simply set to 0，这样当广播操作在扩展维度上运行时，该数组上迭代器的 dataptr 不会移动。

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

数组标量#

另请参阅

标量

数组标量提供了一个 Python 类型层次结构，允许在数组中存储的数据类型和从数组中提取元素时返回的 Python 类型之间建立一对一对应关系。对象数组是一个例外。对象数组是任意 Python 对象的异构集合。当你从对象数组中选择一个项目时，你会得到原始的 Python 对象（而不是对象数组标量，它确实存在，但很少用于实际目的）。

数组标量还提供与数组相同的方法和属性，目的是为了让相同的代码能够用于支持任意维度（包括 0 维）。数组标量是只读的（不可变的），除了 void 标量可以写入，这样结构化数组字段设置可以更自然地工作 (a[0]['f1'] = value)。

索引#

另请参阅

ndarray 上的索引，索引例程

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

  • 整数

  • newaxis

  • 切片

  • Ellipsis

  • 整数数组/类数组（高级）

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

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

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

下一步取决于找到的索引类型。如果所有维度都用整数索引，则返回或设置标量。单个布尔索引数组将调用专门的布尔函数。包含 Ellipsis 或 切片 但不包含高级索引的索引始终会通过计算新的步长和内存偏移量创建对旧数组的视图。然后可以返回此视图，或者对于赋值，可以使用 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\) 个输出。通常在 C 中创建新的通用函数，尽管存在从 Python 函数创建 ufunc 的机制 (frompyfunc)。用户必须提供一个一维循环，该循环实现基本函数，该函数接受输入标量值并将结果标量放置到适当的输出槽中，如实现中所述。

设置#

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

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

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

检查特定于线程的全局变量后，会评估输入以确定 ufunc 应该如何进行，并根据需要构造输入和输出数组。任何不是数组的输入都被转换为数组（如果需要，使用上下文）。哪些输入是标量（因此被转换为 0 维数组）会被记录下来。

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

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

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

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

最后，会做出关于如何执行循环机制的决定，以确保将输入数组的所有元素组合起来，以生成具有正确类型的输出数组。循环执行的选项包括：单循环（对于 :term`连续`、对齐和正确数据类型）、带步长循环（对于非连续但仍对齐且具有正确数据类型）和缓冲循环（对于未对齐或数据类型错误的情况）。根据调用哪个执行方法，设置循环并进行计算。

函数调用#

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

单循环#

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

带步长循环#

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

缓冲循环#

这部分代码处理以下情况：每当输入和/或输出数组未对齐或数据类型错误（包括字节交换）时，底层的一维循环所期望的类型。还假定数组是非连续的。代码的工作原理与带步长循环非常相似，只是内部一维循环进行了修改，以便在 bufsize 块（其中 bufsize 是用户可设置的参数）中对输入执行预处理，对输出执行后处理。底层的一维计算循环在复制过来的数据上被调用（如果需要的话）。在这种情况下，设置代码和循环代码要复杂得多，因为它必须处理

  • 临时缓冲区的内存分配

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

  • 对需要缓冲区的任何输入或输出进行复制和可能的类型转换。

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

  • 将内部 1-D 循环拆分为 bufsize 个块（可能会有余数）。

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

最终输出操作#

通用函数允许其他类似数组的类通过接口无缝传递，因为特定类的输入将导致输出也为该类。这种机制的工作原理如下。如果任何输入不是 ndarray 且定义了 __array_wrap__ 方法，那么具有最大 __array_priority__ 属性的类将决定所有输出的类型（除了传入的任何输出数组）。输入数组的 __array_wrap__ 方法将使用通用函数返回的 ndarray 作为其输入。__array_wrap__ 函数支持两种调用方式。第一种将 ndarray 作为第一个参数，将“上下文”元组作为第二个参数。这是第一个尝试的调用。如果发生 TypeError，则该函数仅使用 ndarray 作为第一个参数调用。

方法#

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

设置#

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

Reduce#

所有通用函数方法都使用相同的底层 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]\) 是 \(i\) 沿着所选轴的第 \(k^{\textrm{th}}\) 个元素。此基本操作会针对具有多个维度的数组重复进行，以便缩减针对沿着所选轴的每个 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*}

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

Reduceat#

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

错位或与输入和/或输出数据类型不匹配的循环数据类型使用缓冲代码处理，其中数据将被复制到临时缓冲区并根据需要转换为正确的数据类型，然后调用底层 1-D 函数。临时缓冲区以不大于用户可设置的缓冲区大小值的（元素）大小创建。因此，循环必须足够灵活，以调用底层 1-D 计算循环足够多次，以便在不大于缓冲区大小的块中完成总计算。

本页内容