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 来移动到下一个元素。下一个元素始终是 C 连续顺序的。该宏通过首先特殊处理 C 连续、1D 和 2D 情况来工作，这些情况非常简单。

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

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

广播#

另请参阅

广播

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

函数 PyArray_Broadcast 采用已定义的迭代器，并使用它们来确定每个维度中的广播形状（要同时创建迭代器并进行广播，请使用函数 PyArray_MultiIterNew）。然后，调整迭代器，使每个迭代器认为它正在迭代具有广播大小的数组。这是通过调整迭代器的维度数以及每个维度中的形状来完成的。这是有效的，因为迭代器步幅也会被调整。广播只调整（或添加）长度为 1 的维度。对于这些维度，步幅变量简单地设置为 0，这样当广播操作遍历扩展的维度时，迭代器的数据指针就不会移动。

在 Numeric 中，广播始终使用扩展维度的零值步幅来实现。在 NumPy 中也是完全一样的方式。主要的区别在于，现在步幅数组被保存在 PyArrayIterObject 中，广播结果涉及的迭代器被保存在 PyArrayMultiIterObject 中，并且调用 PyArray_Broadcast 函数实现了通用广播规则。

数组标量#

另请参阅

标量

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

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

索引#

另请参阅

ndarray 上的索引, 索引例程

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

  • 整数

  • 新轴

  • 切片

  • 省略号

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

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

  • 0-d 布尔值（也包括整数）；0-d 布尔数组是一个特殊情况，必须在高级索引代码中处理。它们表明一个 0-d 布尔数组必须被解释为整数数组。

以及表示整数数组被解释为整数索引的标量数组特殊情况，这很重要，因为整数数组索引强制进行复制，但如果返回标量（完整整数索引）则被忽略。准备好的索引保证有效，除了越界值和高级索引的广播错误。这包括为一个不完整的索引添加一个 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\) 个输出，通过将基本的 1-D 逐元素循环包装成完整易用的函数，这些函数无缝实现广播、类型检查、缓冲转换和输出参数处理。新的通用函数通常用 C 语言创建，尽管有一个机制可以从 Python 函数创建 ufunc（numpy.frompyfunc）。用户必须提供一个 1-D 循环，该循环实现基本函数，接受输入标量值并将结果标量放入适当的输出槽，如实现中所述。

设置#

每次 ufunc 计算都涉及一些与设置计算相关的开销。这种开销的实际意义在于，即使 ufunc 的实际计算速度非常快，你仍然可以编写特定于数组和类型的代码，其速度会比 ufunc 对小数组更快。特别是，使用 ufuncs 对 0-D 数组执行许多计算将比其他基于 Python 的解决方案（默默导入的 scalarmath 模块正是为了让数组标量具有类似 ufunc 的外观和感觉，同时开销显著降低）慢。

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

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

检查完特定于线程的全局变量后，会评估输入以确定 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，以便将其内容复制回输出数组。然后处理输出参数的迭代器。

最后，决定如何执行循环机制，以确保输入数组的所有元素都被组合以生成正确类型的输出数组。循环执行的选项包括单循环（适用于连续、对齐且数据类型正确的数组）、步幅循环（适用于非连续但仍对齐且数据类型正确的数组）以及缓冲循环（适用于错位或数据类型不正确的数组）。根据调用哪种执行方法，然后设置并计算循环。

函数调用#

本节描述了如何为三种不同的执行方式设置和执行基本的通用函数计算循环。如果编译时定义了 NPY_ALLOW_THREADS，并且不涉及对象数组，则在调用循环之前会释放 Python 全局解释器锁 (GIL)。如果需要处理错误条件，则会重新获取它。硬件错误标志仅在 1-D 循环完成后进行检查。

单循环#

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

步幅循环#

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

缓冲循环#

这是处理输入和/或输出数组与底层一维循环期望的数组对齐不当或数据类型错误（包括字节序交换）的情况的代码。数组也被假定为非连续的。这段代码的工作方式非常类似于跨步循环，不同之处在于内部一维循环被修改，以便在 bufsize 块（其中 bufsize 是用户可设置的参数）中对输入进行预处理，并对输出进行后处理。底层一维计算循环被调用到被复制（如果需要）的数据上。在这种情况下，设置代码和循环代码会复杂得多，因为它必须处理

  • 临时缓冲区的内存分配

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

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

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

  • 将内部一维循环分解为 bufsize 块（可能存在余数）。

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

最终输出操作#

通用函数允许其他类数组对象无缝地通过接口传递，因为特定类的输入将促使输出为同一类。实现这一点的机制如下。如果任何输入不是 ndarrays 并且定义了 __array_wrap__ 方法，则具有最大 __array_priority__ 属性的类将决定所有输出的类型（不包括已传入的任何输出数组）。输入数组的 __array_wrap__ 方法将以从通用函数返回的 ndarray 作为其输入而被调用。通用函数支持两种 __array_wrap__ 函数的调用风格。第一种将 ndarray 作为第一个参数，将一个“上下文”元组作为第二个参数。上下文是（ufunc，arguments，output argument number）。这是尝试的第一个调用。如果发生 TypeError，则仅使用 ndarray 作为第一个参数调用该函数。

方法#

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

设置#

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

Reduce#

所有通用函数方法都使用相同的底层一维计算循环，其中输入和输出参数被调整，以便进行适当的约简。例如，reduce 功能的关键在于，一维循环的调用使得输出和第二个输入指向内存中的同一位置，并且两者都具有步长为 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}}\) 个元素。这个基本操作会重复应用于大于一维的数组，以便沿着选定的轴对每个一维子数组进行约简。一个移除选定维度的迭代器处理这个循环。

对于缓冲循环，必须在调用循环函数之前仔细复制和转换数据，因为底层循环期望对齐的数据且数据类型正确（包括字节序）。缓冲循环必须在调用底层函数块（大小不超过用户指定的 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\) 时，每个一维循环处理 \(N\) 个元素。同样，缓冲循环会在调用底层一维计算循环之前复制和转换数据。

Reduceat#

reduceat 函数是 reduce 和 accumulate 函数的泛化。它通过索引指定的输入数组范围来实现 reduce。在循环计算发生之前，会检查额外的索引参数，以确保每个输入在选定维度上不会超出输入数组的大小。循环实现是通过与 reduce 代码非常相似的代码来完成的，该代码重复执行的次数等于索引输入中的元素数量。具体来说：传递给底层一维计算循环的第一个输入指针指向索引数组指示的正确位置的输入数组。此外，传递给底层一维循环的输出指针和第二个输入指针指向内存中的同一位置。一维计算循环的大小固定为当前索引与下一个索引之间的差值（当当前索引是最后一个索引时，下一个索引假定为沿选定维度的数组长度）。这样，一维循环将实现对指定索引的 reduce。

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