NumPy C 代码解释#

狂热主义在于当你忘记目标时加倍努力. — 乔治·桑塔亚纳

权威人士是可以告诉你比你真正想知道的更多关于某事的人. — 未知

该页面试图解释一些新代码背后的逻辑. 这些解释的目的是使人能够更容易地理解实现背后的想法,而不仅仅是盯着代码. 也许通过这种方法,算法可以被更多的人改进,借用和/或优化.

内存模型#

ndarray 的一个基本方面是数组被视为从某个位置开始的内存"块". 对该内存的解释取决于 stride 信息. 对于 \(N\) 维数组中的每个维度,一个整数( stride )指示必须跳过多少字节才能到达该维度中的下一个元素. 除非您有一个单段数组,否则在遍历数组时必须查阅此 stride 信息. 编写接受 strides 的代码并不困难,您只需要使用 char* 指针,因为 strides 的单位是字节. 还要记住,strides 不必是元素大小的单位倍数. 另外,请记住,如果数组的维数为 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 代码中,一个非常常见的操作是需要迭代一个通用的,步长的,N 维数组的所有元素.这种通用 N 维循环的操作被抽象为一个迭代器对象的概念.要编写一个 N 维循环,您只需从一个 ndarray 创建一个迭代器对象,使用迭代器对象结构的 dataptr 成员,并在迭代器对象上调用宏 PyArray_ITER_NEXT 以移动到下一个元素. next 元素始终是 C 连续的顺序.该宏首先特殊处理了 C 连续的,1-D 和 2-D 的情况,这些情况非常简单.

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

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

广播#

参见

广播

在 Numeric (NumPy 的前身) 中,广播的实现在 ufuncobject.c 深处埋藏的几行代码中.在 NumPy 中,广播的概念已经被抽象出来,因此可以在多个地方执行.广播由函数 PyArray_Broadcast 处理.此函数需要传入一个 PyArrayMultiIterObject (或与之二进制等效的东西). PyArrayMultiIterObject 跟踪广播的维度数量和每个维度中的大小,以及广播结果的总大小.它还跟踪正在广播的数组的数量,以及指向正在广播的每个数组的迭代器的指针.

PyArray_Broadcast 函数接收已经定义的迭代器,并使用它们来确定每个维度中的 broadcast 形状(要在广播发生的同时创建迭代器,请使用 PyArray_MultiIterNew 函数).然后,调整迭代器,使每个迭代器都认为它正在迭代具有 broadcast 大小的数组.这是通过调整迭代器的维度数以及每个维度中的 shape 来完成的. 这是可行的,因为迭代器步幅也进行了调整. 广播仅调整(或添加)长度为 1 的维度. 对于这些维度,步幅变量仅设置为 0,这样当广播操作在扩展维度上操作时,该数组上的迭代器的数据指针不会移动.

广播始终使用扩展维度的值为 0 的步幅在 Numeric 中实现. 它在 NumPy 中的实现方式完全相同. 最大的区别在于现在步幅数组保存在 PyArrayIterObject 中,广播结果中涉及的迭代器保存在 PyArrayMultiIterObject 中,并且 PyArray_Broadcast 调用实现了 通用广播规则 .

数组标量#

参见

标量

数组标量提供了一个 Python 类型层次结构,允许数组中存储的数据类型与从数组中提取元素时返回的 Python 类型之间存在一一对应关系.对象数组是这个规则的例外. 对象数组是任意 Python 对象的异构集合. 当您从对象数组中选择一个项目时,您会得到原始的 Python 对象(而不是对象数组标量,它确实存在,但很少用于实际目的).

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

索引#

参见

ndarrays 上的索引 , 索引例程

所有 Python 索引操作 arr[index] 首先通过准备索引并找到索引类型来组织. 支持的索引类型包括:

  • 整数

  • newaxis

  • slice

  • Ellipsis

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

  • 布尔值(单个布尔数组); 如果索引中存在多个布尔数组,或者形状不完全匹配,则布尔数组将转换为整数数组.

  • 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 分配结果.

准备之后,获取和设置相对简单,尽管需要考虑不同的迭代模式.除非在项目获取期间只有一个索引数组,否则会预先检查索引的有效性.否则,为了优化,它在内部循环本身中处理.

通用函数#

参见

通用函数 ( ufunc ) , 通用函数 ( ufunc ) 基础知识

通用函数是可调用对象,它接受 \(N\) 个输入,并通过将逐个元素工作的基本 1-D 循环包装成易于使用的完整函数,从而生成 \(M\) 个输出,这些函数无缝地实现 broadcasting , type-checking , buffered coercion , 和 output-argument handling . 通常在 C 中创建新的通用函数,尽管有一种机制可以从 Python 函数创建 ufuncs( frompyfunc ). 用户必须提供一个 1-D 循环,该循环实现基本函数,该函数获取输入标量值并将结果标量放入适当的输出槽中,如实现中所述.

设置#

每个 ufunc 计算都涉及一些与设置计算相关的开销.这种开销的实际意义在于,即使 ufunc 的实际计算速度非常快,您也可以编写数组和类型特定的代码,这些代码对于小数组来说比 ufunc 工作得更快. 特别是,使用 ufunc 对 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 ,以便将其内容复制回输出数组.然后处理输出参数的迭代器.

最后,决定如何执行循环机制,以确保输入数组的所有元素都组合起来以生成正确类型的输出数组.循环执行的选项包括单循环(对于 contiguous ,对齐和正确的数据类型),跨步循环(对于非连续但仍然对齐和正确的数据类型)和缓冲循环(对于未对齐或不正确的数据类型情况).根据调用哪种执行方法来设置和计算循环.

函数调用#

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

单循环#

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

跨步循环#

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

缓冲循环#

每当输入和/或输出数组未对齐或数据类型错误(包括字节交换)时,此代码都会处理这种情况,这与底层1-D循环所期望的不同. 数组也被假定为非连续的. 除了内部 1-D 循环被修改为对输入执行预处理以及对 bufsize 大小的输出执行后处理之外(其中 bufsize 是用户可设置的参数),此代码的工作方式与跨步循环非常相似. 底层的 1-D 计算循环在复制过来的数据上调用(如果需要). 在这种情况下,设置代码和循环代码要复杂得多,因为它必须处理:

  • 临时缓冲区的内存分配

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

  • 复制和可能转换任何需要缓冲区的输入或输出的数据.

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

  • 将内部 1-D 循环分解为 bufsize 大小的块(可能带有余数).

同样,硬件错误标志在每个 1-D 循环结束时进行检查.

最终输出操作#

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

方法#

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

设置#

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

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]\) 是沿所选轴的 \(k^{\textrm{th}}\) 的第 \(i\) 个元素.对于具有大于 1 维的数组,此基本操作会重复执行,以便对沿所选轴的每个一维子数组执行归约.删除所选维度的迭代器处理此循环.

对于缓冲循环,必须在调用循环函数之前复制和强制转换数据,因为底层循环需要正确数据类型的对齐数据(包括字节顺序).缓冲循环必须在 bufsize 不超过用户指定的 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\) 时,每个一维循环对 \(N+1\) 个元素进行操作.同样,缓冲循环会注意在调用底层一维计算循环之前复制和强制转换数据.

Reduceat#

reduceat 函数是 reduce 和 accumulate 函数的概括.它实现了由索引指定的输入数组范围上的 reduce .在循环计算发生之前,会检查额外的 indices 参数,以确保每个输入对于沿所选维度的输入数组而言不会太大.循环实现使用与 reduce 代码非常相似的代码进行处理,该代码根据 indices 输入中的元素数量重复多次.特别是:传递到基础一维计算循环的第一个输入指针指向索引数组指示的正确位置处的输入数组.此外,传递到基础一维循环的输出指针和第二个输入指针指向内存中的相同位置.一维计算循环的大小固定为当前索引和下一个索引之间的差(当当前索引是最后一个索引时,则假定下一个索引是沿所选维度的数组的长度).这样,一维循环将实现指定索引上的 reduce .

未对齐或与输入和/或输出数据类型不匹配的循环数据类型使用缓冲代码处理,其中数据被复制到临时缓冲区,并在调用底层一维函数之前根据需要强制转换为正确的数据类型.临时缓冲区以不大于用户可设置的 buffer-size 值的(元素)大小创建.因此,循环必须足够灵活,才能多次调用底层一维计算循环,以完成不大于缓冲区大小的块中的总计算.

On this page