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-D 和 2-D 的情况，这些情况的工作方式非常简单。

对于一般情况，迭代的工作方式是在迭代器对象中跟踪坐标计数器列表。在每次迭代时，最后一个坐标计数器都会增加（从 0 开始）。如果此计数器小于该维度中数组大小减 1（预先计算和存储的值），则增加计数器，并且 dataptr 成员按该维度的步长增加，宏结束。如果到达维度的末尾，则最后一个维度的计数器重置为零，并且通过减去步长值乘以该维度中元素数量减 1，将 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 中，广播始终使用扩展维度的 0 值步长来实现。在 NumPy 中，它以完全相同的方式完成。最大的区别在于，现在步长数组被保存在 PyArrayIterObject 中，广播结果中涉及的迭代器被保存在 PyArrayMultiIterObject 中，而 PyArray_Broadcast 调用实现了通用广播规则。

数组标量#

另请参阅

标量

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

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

索引#

另请参阅

ndarray 的索引, 索引例程

所有 Python 索引操作 arr[index] 的组织方式是首先准备索引并查找索引类型。支持的索引类型为

  • 整数

  • newaxis

  • 切片

  • 省略号

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

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

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

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

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

高级索引#

到目前为止，最复杂的情况是高级索引，它可能会也可能不会与典型的基于视图的索引相结合。在尝试理解这一点之前，您可能需要熟悉它的细微之处。高级索引代码有三个不同的分支和一个特殊情况

  • 有一个索引数组，并且它以及赋值数组可以轻松地迭代。例如，它们可能是连续的。此外，索引数组的类型必须为 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 对象，但不会被初始化为对象，因为它仅在内部使用）。此循环对象具有与 PyArray_Broadcast 一起使用的布局，以便广播可以像在代码的其他部分中一样处理。

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

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

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

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

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

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

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

函数调用#

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

单循环#

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

跨步循环#

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

缓冲循环#

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

  • 临时缓冲区的内存分配

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

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

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

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

同样，硬件错误标志会在每个一维循环的末尾进行检查。

最终输出操作#

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

方法#

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

设置#

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

Reduce#

所有 ufunc 方法都使用相同的底层一维计算循环，并调整输入和输出参数，以便进行适当的归约。例如，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}}\) 个元素。对于大于 1 维的数组重复此基本操作，以便在沿所选轴的每个一维子数组上进行归约。删除所选维度的迭代器处理此循环。

对于缓冲循环，必须注意在调用循环函数之前复制和转换数据，因为底层循环需要正确数据类型（包括字节顺序）的对齐数据。缓冲循环必须在以不大于用户指定的 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。

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