NumPy C 代码解释#
狂热主义在于当你忘记了自己的目标时,加倍努力. — 乔治·桑塔亚纳
权威是指一个人,他可以告诉你一些你其实不太想知道的事情. — 未知
本页面试图解释一些新代码背后的逻辑.这些解释的目的是使某人能够比仅仅盯着代码更容易理解实现背后的思想.也许通过这种方式,算法可以得到改进、借鉴和/或优化的人会更多.
内存模型#
一个基本方面是 ndarray 被视为从某个位置开始的”内存块”.对这段内存的解释取决于 步幅 信息.在 \(N\) 维数组中的每个维度,一个整数（步幅）规定了必须跳过多少字节才能到达该维度中的下一个元素.除非你有一个单段数组,否则在遍历数组时必须查阅这个 步幅 信息.编写接受步幅的代码并不难,你只需要使用 char* 指针,因为步幅是以字节为单位的.还要记住,步幅不必是元素大小的单位倍数.此外,记住如果数组的维数为 0（有时称为 rank-0 数组）,那么 步幅 和 维度 变量是 NULL.
除了包含在 PyArrayObject 的 strides 和 dimensions 成员中的结构信息外,flags 还包含关于如何访问数据的重要信息.特别是,当内存位于适合数据类型数组的边界时,会设置 NPY_ARRAY_ALIGNED 标志.即使你有一个 连续 的内存块,也不能假设直接解引用指向元素的数据类型指针是安全的.只有在设置了 NPY_ARRAY_ALIGNED 标志时,这才是安全的操作.在某些平台上它可能会工作,但在其他平台上,如 Solaris,会导致总线错误.如果你计划写入数组的内存区域,还应确保 NPY_ARRAY_WRITEABLE 标志.也可以获取指向不可写内存区域的指针.有时,当未设置 NPY_ARRAY_WRITEABLE 标志时写入内存区域只是不礼貌.其他时候可能会导致程序崩溃（例如,一个只读的内存映射文件的数据区域）.
数据类型封装#
The 数据类型 是 ndarray 的重要抽象.操作将查看数据类型以提供操作数组所需的关键功能.此功能由 PyArray_Descr 结构体的 f 成员指向的函数指针列表提供.通过这种方式,只需在 f 成员中提供带有合适函数指针的 PyArray_Descr 结构体,即可扩展数据类型的数量.对于内置类型,有一些优化会绕过此机制,但数据类型抽象的目的是允许添加新的数据类型.
内置数据类型之一,:class:void 数据类型允许包含1个或多个字段的任意 结构化数据类型 作为数组的元素.:term:字段 仅仅是另一个数据类型对象,以及当前结构化类型中的偏移量.为了支持任意嵌套的字段,为void类型实现了几种递归的数据类型访问实现.一个常见的习惯用法是遍历字典的元素,并根据存储在给定偏移量的数据类型对象执行特定操作.这些偏移量可以是任意数字.因此,必须认识到并考虑到遇到未对齐数据的可能性（如果必要的话）.
N-D 迭代器#
参见
在许多NumPy代码中非常常见的操作是需要遍历一个通用的、跨步的、N维数组的所有元素.这种通用N维循环的操作在迭代器对象的概念中被抽象化.要编写一个N维循环,你只需要从一个ndarray创建一个迭代器对象,使用迭代器对象结构的 dataptr 成员,并在迭代器对象上调用宏 PyArray_ITER_NEXT 以移动到下一个元素.``next`` 元素总是按C连续顺序排列.该宏首先通过特殊处理C连续、1-D和2-D情况来工作,这些情况非常简单.
对于一般情况,迭代通过在迭代器对象中跟踪一个坐标计数器列表来工作.在每次迭代中,最后一个坐标计数器会增加（从0开始）.如果这个计数器小于该维度数组大小减一（这是一个预先计算并存储的值）,那么计数器会增加,并且 dataptr 成员会按该维度的步幅增加,宏结束.如果到达一个维度的末尾,最后一个维度的计数器重置为零,并且 dataptr 通过减去步幅值乘以该维度元素数减一（这也是预先计算并存储在 backstrides 成员中的值）回到该维度的开始.在这种情况下,宏不会结束,但一个局部维度计数器会递减,以便倒数第二个维度取代最后一个维度的角色,并在倒数第二个维度上再次执行前面描述的测试.通过这种方式,:c:member:dataptr 会根据任意步幅适当调整.
在 PyArrayIterObject 结构体中,除非底层数组是 C 连续的,否则 coordinates 成员会维护当前的 N 维计数器,在这种情况下,坐标计数会被绕过.:c:member:index 成员会跟踪迭代器的当前扁平索引.它由 PyArray_ITER_NEXT 宏更新.
广播#
参见
在 Numeric,NumPy 的前身,广播是在 ufuncobject.c 深处几行代码中实现的.在 NumPy 中,广播的概念已经被抽象化,以便可以在多个地方执行.广播由函数 PyArray_Broadcast 处理.这个函数需要一个 PyArrayMultiIterObject`（或等效的二进制对象）作为输入.:c:type:`PyArrayMultiIterObject 跟踪广播的维度数量和每个维度的大小,以及广播结果的总大小.它还跟踪正在广播的数组数量,并为每个正在广播的数组维护一个指向迭代器的指针.
The PyArray_Broadcast function takes the iterators that have already
been defined and uses them to determine the broadcast shape in each
dimension (to create the iterators at the same time that broadcasting
occurs then use the PyArray_MultiIterNew function).
Then, the iterators are
adjusted so that each iterator thinks it is iterating over an array
with the broadcast size. This is done by adjusting the iterators
number of dimensions, and the shape in each dimension. This works
because the iterator strides are also adjusted. Broadcasting only
adjusts (or adds) length-1 dimensions. For these dimensions, the
strides variable is simply set to 0 so that the data-pointer for the
iterator over that array doesn’t move as the broadcasting operation
operates over the extended dimension.
在 Numeric 中,广播总是通过扩展维度使用 0 值的步幅来实现的.在 NumPy 中也是以完全相同的方式完成的.最大的区别在于,现在步幅数组在 PyArrayIterObject 中被跟踪,广播结果中涉及的迭代器在 PyArrayMultiIterObject 中被跟踪,并且 PyArray_Broadcast 调用实现了 通用广播规则.
数组标量#
参见
数组.标量
数组标量提供了一个层次结构的 Python 类型,允许存储在数组中的数据类型与从数组中提取元素时返回的 Python 类型之间的一对一对应关系.这个规则的一个例外是对象数组.对象数组是任意 Python 对象的异构集合.当你从对象数组中选择一个项目时,你会得到原始的 Python 对象（而不是存在的对象数组标量,但在实际用途上很少使用）.
数组标量也提供了与数组相同的方法和属性,目的是相同的代码可以用来支持任意维度（包括0维）.数组标量是只读的（不可变的）,除了void标量也可以写入,以便结构化数组字段设置更自然（a[0]['f1'] = value）.
索引#
参见
所有 Python 索引操作 arr[index] 都是通过首先准备索引并找到索引类型来组织的.支持的索引类型有:
整数
整数数组/类数组对象 (高级)
布尔值（单个布尔数组）;如果有多个布尔数组作为索引或形状不完全匹配,布尔数组将转换为整数数组.
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 被要求分配结果,否则必须由调用者处理.
在准备之后,获取和设置相对简单,尽管需要考虑不同的迭代模式.除非在获取项目时只有一个索引数组,否则会事先检查索引的有效性.否则,为了优化,会在内部循环中处理.
通用函数#
通用函数是可调用对象,它们通过将基本的1-D循环（逐元素工作）封装成易于使用的完整函数,无缝实现 广播、类型检查、缓冲强制 和 输出参数处理,从而接受 \(N\) 个输入并产生 \(M\) 个输出.新的通用函数通常在C中创建,尽管有一种机制可以从Python函数创建ufuncs（frompyfunc）.用户必须提供一个1-D循环,该循环实现基本函数,接受输入标量值并将结果标量放入适当的输出槽中,如实现中所述.
设置#
每个 ufunc 计算都涉及一些与设置计算相关的开销.这种开销的实际意义在于,即使 ufunc 的实际计算非常快,你仍然能够编写针对数组和类型特定的代码,这些代码对于小数组来说会比 ufunc 更快.特别是,使用 ufuncs 对许多 0-D 数组进行计算会比其他基于 Python 的解决方案更慢（静默导入的 scalarmath 模块的存在正是为了使数组标量具有与 ufunc 计算相似的外观和感觉,同时显著减少开销）.
当调用 ufunc 时,必须完成许多事情.从这些设置操作中收集的信息存储在一个循环对象中.这个循环对象是一个 C 结构（可能会变成一个 Python 对象,但不会初始化为这样的对象,因为它仅在内部使用）.这个循环对象具有与 PyArray_Broadcast 一起使用的布局,以便广播可以以与其他代码部分相同的方式处理.
首先做的事情是在线程特定的全局字典中查找当前的缓冲区大小、错误掩码和关联的错误对象的值.错误掩码的状态控制当发现错误条件时会发生什么.需要注意的是,只有在每次执行1-D循环后才会检查硬件错误标志.这意味着如果输入和输出数组是连续的并且类型正确,以至于执行单个1-D循环,那么标志可能直到数组的所有元素都被计算后才会被检查.在线程特定的字典中查找这些值需要时间,除了非常小的数组外,这很容易被忽略.
在检查线程特定的全局变量后,输入被评估以确定ufunc应如何进行,并且在必要时构建输入和输出数组.任何不是数组的输入都被转换为数组（如有必要使用上下文）.注意哪些输入是标量（因此被转换为0-D数组）.
接下来,根据输入数组类型从可用的 1-D 循环中选择一个适当的 1-D 循环用于 ufunc.这个 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 以便将其内容复制回输出数组.随后处理输出参数的迭代器.
最后,关于如何执行循环机制的决定被做出,以确保输入数组的所有元素都被组合以生成正确类型的输出数组.循环执行的选项包括单循环（用于 :term`连续` 、对齐且数据类型正确的情况）、跨步循环（用于非连续但仍对齐且数据类型正确的情况）和缓冲循环（用于未对齐或数据类型不正确的情况）.根据调用的执行方法,循环随后被设置并计算.
函数调用#
本节描述了如何为三种不同类型的执行设置和执行基本的通用函数计算循环.如果在编译期间定义了 NPY_ALLOW_THREADS ,那么在不涉及对象数组的情况下,Python 全局解释器锁 (GIL) 会在调用循环之前释放.如果需要处理错误条件,则会重新获取 GIL.硬件错误标志仅在完成 1-D 循环后检查.
一个循环#
这是最简单的情况.ufunc 通过调用底层 1-D 循环一次来执行.只有在输入和输出的数据类型（包括字节顺序）对齐且所有数组具有均匀步长（连续、0-D 或 1-D）时,才可能实现这一点.在这种情况下,1-D 计算循环被调用一次来计算整个数组的计算.请注意,硬件错误标志仅在整个计算完成后才会被检查.
Strided loop#
当输入和输出数组对齐且类型正确,但步幅不均匀（非连续且为2-D或更大）时,则使用第二个循环结构进行计算.这种方法将所有输入和输出参数的迭代器转换为遍历除最大维度外的所有维度.内部循环由底层的1-D计算循环处理.外部循环是转换后的迭代器上的标准迭代器循环.每次完成1-D循环后都会检查硬件错误标志.
Buffered loop#
这是处理输入和/或输出数组与底层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, 参数, 输出参数编号).这是首先尝试的调用.如果发生 TypeError,则该函数仅以 ndarray 作为第一个参数进行调用.
方法#
有三种 ufuncs 方法需要类似于通用 ufuncs 的计算.这些是 ufunc.reduce、ufunc.accumulate 和 ufunc.reduceat.每种方法都需要一个设置命令,然后是一个循环.对于这些方法,有四种可能的循环样式,分别对应于无元素、单元素、跨步循环和缓冲循环.这些是与通用函数调用实现的基本循环样式相同,除了无元素和单元素情况,这些是特殊情况,分别发生在输入数组对象有 0 和 1 个元素时.
设置#
所有三种方法的设置函数是 construct_reduce.这个函数创建一个减少循环对象,并用完成循环所需的参数填充它.所有方法仅适用于接受2个输入并返回1个输出的ufuncs.因此,假设签名是 [otype, otype, otype],其中 otype 是请求的减少数据类型,选择底层1-D循环.然后从（每个线程）全局存储中检索缓冲区大小和错误处理.对于未对齐或数据类型不正确的小数组,会进行复制,以便使用未缓冲的部分代码.然后,选择循环策略.如果数组中有1个或0个元素,则选择简单的循环方法.如果数组未对齐且具有正确的数据类型,则选择跨步循环.否则,必须执行缓冲循环.然后,建立循环参数,并构造返回数组.输出数组的 形状 取决于方法是否是 reduce、accumulate 或 reduceat.如果已经提供了输出数组,则检查其形状.如果输出数组不是C连续的、对齐的和正确的数据类型,则使用设置了 NPY_ARRAY_WRITEBACKIFCOPY 标志的临时副本.这样,方法将能够使用行为良好的输出数组,但在函数完成时调用 PyArray_ResolveWritebackIfCopy 时,结果将被复制回真正的输出数组.最后,设置迭代器以循环正确的 :term:`轴`（取决于提供给方法的轴值）,并且设置例程返回到实际计算例程.
Reduce#
所有ufunc方法使用相同的底层1-D计算循环,输入和输出参数经过调整,以便进行适当的归约.例如,:meth:reduce 功能的关键在于,1-D循环在调用时,输出和第二个输入指向内存中的同一位置,并且两者都有0的步长.第一个输入指向输入数组,步长由所选轴的适当步幅给出.通过这种方式,执行的操作是
其中 \(N+1\) 是输入中的元素数量,:math:i,:math:o 是输出,而 \(i[k]\) 是沿所选轴的 \(i\) 的第 \(k^{\textrm{th}}\) 个元素.对于维度大于1的数组,此基本操作会重复进行,以便沿所选轴对每个1-D子数组进行归约.移除了所选维度的迭代器处理此循环.
对于缓冲循环,必须在调用循环函数之前复制和转换数据,因为底层循环期望对齐的正确数据类型数据（包括字节顺序）.缓冲循环必须在调用循环函数之前处理这种复制和转换,处理的数据块大小不超过用户指定的 bufsize.
累积#
accumulate 方法与 reduce 方法非常相似,因为输出和第二个输入都指向输出.不同之处在于,第二个输入指向当前输出指针后面一个步长的内存.因此,执行的操作是
输出与输入具有相同的形状,并且每个1-D循环在选定轴的形状为 \(N+1\) 时操作 \(N\) 个元素.同样,缓冲循环在调用底层的1-D计算循环之前会复制和转换数据.
Reduceat#
reduceat 函数是 reduce 和 accumulate 函数的泛化.它在输入数组指定索引范围内的实现 reduce .额外的索引参数被检查以确保每个输入对于输入数组在选定维度上不会太大,然后进行循环计算.循环实现使用与 reduce 代码非常相似的代码处理,重复次数与索引输入中的元素数量相同.特别是:传递给底层 1-D 计算循环的第一个输入指针指向输入数组在索引数组指示的正确位置.此外,传递给底层 1-D 循环的输出指针和第二个输入指针指向内存中的同一位置.1-D 计算循环的大小固定为当前索引和下一个索引之间的差值（当当前索引是最后一个索引时,下一个索引假定为数组在选定维度上的长度）.通过这种方式,1-D 循环将在指定索引上实现 reduce .
如果对齐不正确或循环数据类型与输入和/或输出数据类型不匹配,则使用缓冲代码处理,其中数据在调用底层一维函数之前被复制到临时缓冲区,并在必要时转换为正确的数据类型.临时缓冲区以（元素）大小创建,不大于用户可设置的缓冲区大小值.因此,循环必须足够灵活,能够调用底层一维计算循环足够多次,以分块完成总计算,且每个块不大于缓冲区大小.