NumPy 核心 C 实现原理深度解析内存模型、迭代器、广播与 ufunc 全流程【免费下载链接】numpyThe fundamental package for scientific computing with Python.项目地址: https://gitcode.com/gh_mirrors/nu/numpy本篇指南以 NumPy 官方开发文档 internals.code-explanations.rst 为主体系统讲解 ndarray 的内存布局stride 与 flags、数据类型封装、N 维迭代器、广播机制、数组标量、索引含高级索引以及通用函数ufunc从设置到执行的全流程。读完本文你将理解 NumPy C 层的关键设计思想——包括 strided 内存访问为何必须用char*、广播为什么用零步长实现、ufunc 的三种循环执行策略如何选择——并能沿着源码路径深入阅读 numpy/_core/src/multiarray/ 与 numpy/_core/src/umath/ 中的真实实现。本文的目的正如原文档开篇所言解释新代码背后的逻辑让读者比死盯代码更容易理解实现思路从而能够改进、借鉴并优化这些算法。以下各节均以官方文档为骨架并补充了当前仓库中可核对的源码位置与实现细节。内存模型Memory model:class:ndarray最根本的特征是数组被看作从某个起始位置开始的一段内存块。这段内存如何解释完全取决于stride步长信息。对每一个维度数组都用一个整数stride表示要跳过多少字节才能到达该维度的下一个元素。对于 N 维数组遍历时需要依据 stride 信息移动指针除非数组是单段single-segment连续数组。关键要点编写接受 stride 的代码必须使用char*指针因为 stride 的单位是字节stride 不必是元素大小的整数倍例如广播或某些切片视图会出现非整数倍的步长如果数组维度数为 0即所谓的rank-0数组则 strides 和 dimensions 变量为NULL。除了保存在PyArrayObject的 strides 与 dimensions 成员中的结构信息外flags 标志位还包含数据访问方式的重要信息NPY_ARRAY_ALIGNED当内存位于符合 dtype 要求的合适边界上时置位。即使你拥有一段连续内存也不能想当然地假设按数据类型解引用指针是安全的——只有当NPY_ARRAY_ALIGNED置位时才是安全操作。在某些平台上碰巧能工作但在另一些平台如 Solaris上会触发总线错误bus error。该标志位在源码中广泛参与判定例如 numpy/_core/src/multiarray/array_coercion.c 在数组创建/转换路径中维护这些标志。NPY_ARRAY_WRITEABLE如果计划写入数组内存区域必须确保该标志置位。有时拿到的是不可写内存区域的指针——写入时可能只是不礼貌有时则会导致程序崩溃例如数据区域是只读的内存映射文件。从源码结构看这些标志位贯穿了数组创建ctors.c、转换convert.c和 ufunc 调度ufunc_object.c等几乎所有路径。数据类型封装Data-type encapsulation:ref:datatype arrays.dtypes是 ndarray 的重要抽象层操作会依赖 datatype 来提供操作数组所需的关键功能。这些功能通过PyArray_Descr结构中f成员指向的一串函数指针提供。由此带来的扩展性至关重要只要提供一个带有合适函数指针位于f成员的PyArray_Descr结构就能新增一种数据类型。内置类型有一些绕过该机制的优化但该抽象的根本目的就是允许添加新类型。关于 void结构化类型的递归处理内置数据类型之一——void 类型——允许数组元素是包含 1 个或多个字段的任意 :term:structured types structured data type。字段field本质上是另一个 datatype 对象 一个相对当前结构化类型的偏移量。为支持任意嵌套字段void 类型实现了多个递归的数据类型访问实现。常见惯用法是遍历字典中的元素根据给定偏移处存储的 datatype 对象执行特定操作。这些偏移量可以是任意数值因此必须意识到可能遇到未对齐数据并在必要时加以处理。更完整的用户视角说明可参考 doc/source/reference/arrays.dtypes.rst。N 维迭代器N-D iterators遍历一个通用、带 stride 的 N 维数组的全部元素是 NumPy 代码中极为常见的操作。这一通用 N 维循环被抽象为迭代器对象只需从 ndarray 创建迭代器对象使用迭代器结构中的dataptr成员调用宏PyArray_ITER_NEXT移动到下一个元素——下一个元素永远按 C 连续C-contiguous顺序。宏的实现策略首先对 C 连续、1 维和 2 维这些简单情况做特判因为这些情况非常简单。一般情况下迭代通过维护迭代器对象中的坐标计数器列表coordinate counters完成每次迭代最后一个坐标计数器递增从 0 开始若该计数器小于该维数组大小减 1预计算并存储的值则递增计数器dataptr按该维 stride 增加宏结束若到达维度末尾最后一个维度的计数器重置为 0dataptr通过减去strides * (元素数 - 1)也预计算并存储在迭代器对象的backstrides成员中移回该维起点此时宏不结束而是递减一个局部维度计数器让倒数第二维接替最后一维的角色重复上述测试——由此dataptr可以适应任意 striding。另外两个成员coordinates维护当前 N 维计数器除非底层数组是 C 连续的此时坐标计数被绕过index跟踪当前扁平索引由PyArray_ITER_NEXT宏更新。在仓库中backstrides的预计算可见于 numpy/_core/src/multiarray/iterators.c例如it-backstrides[i] it-strides[i] * it-dims_m1[i]PyArray_ITER_NEXT被 numpy/_core/src/umath/ufunc_object.c、numpy/_core/src/multiarray/mapping.c 等核心路径大量调用。用户侧对应的就是numpy.nditer参见 doc/source/reference/arrays.nditer.rst。广播Broadcasting在 NumPy 的前身 Numeric 中广播只实现为深埋在ufuncobject.c里的几行代码。而在 NumPy 中广播被抽象出来可以在多个地方执行。核心函数是PyArray_Broadcast其实现位于 numpy/_core/src/multiarray/iterators.c#L1147。它要求传入一个PyArrayMultiIterObject或其二进制等价物。该对象跟踪广播后的维度数nd与每个维度的大小dimensions广播结果的总大小size参与广播的数组个数numiter每个参与广播数组的迭代器指针数组iters。PyArray_Broadcast的工作流程与源码逐行对应发现广播维度数取所有参与数组ndim的最大值nd PyArray_MAX(...)发现每一维的广播形状遍历每个维度先把dimensions[i]初始化为 1再对每个数组取对应维的大小tmp若tmp 1则跳过长度 1 的维度可被广播否则若当前值为 1 则采纳tmp若已非 1 且不等于tmp则抛出 shape mismatch 异常set_shape_mismatch_exception重置每个迭代器的维度与步长使用 0 值步长实现广播源码注释 using 0 valued strides for broadcasting。对新增的维度或底层数组该维大小为 1 的维度将it-strides[j] 0从而在广播操作沿扩展维度进行时该数组的数据指针不动否则沿用原数组步长并更新backstrides与factors。广播只调整或增加长度为 1 的维度——这些维度的 stride 直接设为 0。这与 Numeric 的做法完全一致Numeric 一直用 0 值 stride 表示扩展维度区别在于现在步长数组保存在PyArrayIterObject中、参与广播的迭代器由PyArrayMultiIterObject跟踪PyArray_Broadcast调用实现了通用的广播规则。若想在创建迭代器的同时完成广播则使用PyArray_MultiIterNew函数。用户侧规则说明参见 doc/source/user/basics.broadcasting.rst。数组标量Array scalars数组标量提供了一组 Python 类型层级使数组中存储的数据类型与从数组提取元素时返回的 Python 类型形成一一对应关系。例外是对象数组对象数组是任意 Python 对象的异构集合从中取出的元素返回原始 Python 对象而不是对象数组标量——虽然它确实存在但很少实际使用。数组标量还提供与数组相同的方法和属性意图是让同一套代码可以支持任意维度包括 0 维。数组标量是**只读不可变**的唯二例外是 void 标量——它也可以被写入以便结构化数组字段赋值a[0][f1] value更自然。详细说明参见 doc/source/reference/arrays.scalars.rst。索引Indexing所有 Python 索引操作arr[index]的组织方式是先准备索引并判定索引类型。支持的索引类型包括整数integernewaxissliceEllipsis整数数组/类数组高级索引advanced布尔数组单个布尔数组若索引中有多个布尔数组或形状不完全匹配布尔数组会被转换为整数数组0 维布尔以及整数数组0 维布尔数组是高级索引代码中必须特殊处理的情形它们标记0 维布尔数组必须被解释为整数数组此外还有标量数组特例标记整数数组被解释为整数索引。这很重要因为整数数组索引强制产生副本但若返回标量完整整数索引则被忽略。准备好的索引保证有效唯一例外是越界值和高级索引的广播错误。这包括对不完整的索引会自动补上Ellipsis例如用单个整数索引二维数组时。下一步取决于找到的索引类型若所有维度都被整数索引则返回或设置一个标量单个布尔索引数组会调用专门的布尔函数包含Ellipsis或slice但没有高级索引的索引总是通过计算新步长和内存偏移在原数组上创建视图该视图可直接返回或赋值时用PyArray_CopyObject填充。注意在对象dtype数组的复杂赋值中PyArray_CopyObject也可能在其他分支的临时数组上被调用其实现位于 numpy/_core/src/multiarray/array_coercion.c 等文件中。高级索引Advanced indexing高级索引是最复杂的情形它可能与传统基于视图的索引结合。这里整数索引被解释为基于视图的。高级索引代码有三个分支和一个特例只有一个索引数组且它以及赋值数组可以平凡地迭代例如它们是连续的且索引数组必须是intp类型赋值时的值数组类型要正确。这纯粹是快速路径只有整数数组索引不存在子数组视图索引与高级索引混合视图索引定义了一组子数组由高级索引组合。例如arr[[1, 2, 3], :]通过垂直堆叠子数组arr[1, :]、arr[2, :]、arr[3, :]创建存在子数组但恰好只有一个元素这种情况可按无子数组处理但设置阶段需要小心。判定属于哪种情形、检查广播、确定所需的转置类型全部在PyArray_MapIterNew中完成源码位于 numpy/_core/src/multiarray/mapping.c 及其头文件 numpy/_core/src/multiarray/mapping.h。设置之后有两种情况无子数组或子数组只有一个元素无需子数组迭代准备一个迭代器同时迭代所有索引数组以及结果/值数组存在子数组准备三个迭代器——一个用于索引数组、一个用于结果/值数组减去其子数组部分、一个用于原数组与结果/赋值数组的子数组。前两个迭代器给出或允许计算出子数组起始位置的指针从而可以重启子数组迭代。当高级索引彼此相邻时可能需要转置。所有必要的转置由PyArray_MapIterSwapAxes处理除非PyArray_MapIterNew被要求分配结果否则调用方必须处理转置。准备完成后获取get与设置set相对直接但需要考虑不同的迭代模式。除非在取值时只有一个索引数组否则索引的有效性会事先检查否则为了优化在内部循环本身处理。用户侧文档参见 doc/source/reference/arrays.indexing.rst 与 doc/source/user/basics.indexing.rst。通用函数Universal functions / ufunc通用函数ufunc是可调用对象它们接收 N 个输入、产生 M 个输出将逐元素工作的基本 1-D 循环包装成易用的完整函数并无缝实现广播、类型检查casting、缓冲强制转换和输出参数处理。新 ufunc 通常在 C 中创建但也有从 Python 函数创建 ufunc 的机制frompyfunc。用户必须提供一个实现基本功能的 1-D 循环将输入标量值计算结果放入相应输出槽位。完整说明参见 doc/source/reference/ufuncs.rst。设置Setup每次 ufunc 计算都有设置开销。实际意义是即使 ufunc 的实际计算非常快针对小数组手写数组/类型特定的代码会比 ufunc 更快。特别是用 ufunc 对 0 维数组做大量计算会比纯 Python 方案慢——静默导入的scalarmath模块正是为了以显著更低的开销给数组标量提供 ufunc 计算的外观与手感。调用 ufunc 时需要完成多件事收集到的信息存放在一个**循环对象loop object**中——这是一个 C 结构体本可成为 Python 对象但因为仅内部使用而未被如此初始化。该对象的布局兼容PyArray_Broadcast因此广播可以用与其他代码段相同的方式处理。设置流程对照源码 numpy/_core/src/umath/ufunc_object.c 中PyUFunc_GenericFunctionInternal的实现读取线程局部全局字典获取当前的 buffer-size缓冲大小、error mask错误掩码和关联的错误对象。错误掩码的状态决定发现错误条件时的行为。注意硬件错误标志只在每个 1-D 循环执行之后检查——这意味着若输入输出数组连续且类型正确只执行单个 1-D 循环时标志可能要等到整个数组算完才被检查。该步在源码中由_get_bufsize_errmask完成同时取出buffersize与errormask。查找线程特定字典需要时间但对除极小数组外的所有情况都可忽略检查线程局部全局变量后评估输入决定 ufunc 如何继续必要时构造输入/输出数组。所有非数组输入被转换为数组必要时使用 context并记录哪些输入是标量因此被转换为 0 维数组选择 1-D 循环根据输入数组类型从 ufunc 可用的 1-D 循环中选择——尝试把输入的数据类型签名与可用签名匹配。内置类型的签名存储在 ufunc 结构的ufunc.types成员中用户自定义类型的签名存储在以参数列表中第一个用户自定义类型的类型号为键的userloops字典头部元素为CObject所指向的函数信息链表中。签名搜索持续到找到所有输入数组都能安全转换的签名忽略不允许决定结果类型的标量参数。搜索过程的含义存储签名时应把较小类型放在较大类型之下。若找不到 1-D 循环则报错否则用存储的签名更新argument_list——以防需要转换并固定 1-D 循环假定的输出类型特殊检查若 ufunc 有 2 个输入、1 个输出且第二个输入是Object数组则执行特殊检查当第二个输入不是 ndarray、拥有__array_priority__属性且有__r{op}__特殊方法时返回NotImplemented。这样 Python 会给予另一个对象完成操作的机会而不是使用通用的对象数组计算——例如这允许稀疏矩阵覆盖乘法运算符的 1-D 循环输入处理对于小于指定缓冲大小的输入数组对所有非连续、未对齐或字节序不对的数组制作副本确保小数组使用单个循环然后为所有输入数组创建数组迭代器并将迭代器集合广播到单一形状输出处理处理输出参数若有构造缺失的返回数组。若提供的输出数组类型不对或未对齐且小于缓冲大小则构造带NPY_ARRAY_WRITEBACKIFCOPY标志的新输出数组函数结束时调用PyArray_ResolveWritebackIfCopy将其内容复制回输出数组。然后处理输出参数的迭代器选择循环执行机制决定如何执行循环以组合所有输入元素并产生正确类型的输出。三种选择one-loop连续、对齐且类型正确的数据、strided-loop非连续但仍对齐且类型正确、buffered loop未对齐或类型错误的情形。依据所选执行方法设置并计算循环。函数调用Function call本小节描述三种执行方式下基本通用函数计算循环的设置与执行。若编译时定义了NPY_ALLOW_THREADS则只要不涉及对象数组调用循环前会释放 Python GIL全局解释器锁处理错误条件需要时再重新获取。硬件错误标志只在 1-D 循环完成后检查。单循环One loop最简单的情形ufunc 通过恰好调用一次底层 1-D 循环执行。只有输入输出都是对齐的、类型含字节序正确的数据且所有数组步长均匀连续、0-D 或 1-D时才可能。此时硬件错误标志在整个计算完成后才检查。跨步循环Strided loop当输入输出数组对齐且类型正确但步长不均匀非连续且 2 维或更高时使用第二种循环结构把输入/输出参数的所有迭代器转换为除最大维度外全部迭代的形式内层循环交给底层 1-D 计算循环外层循环是转换后迭代器上的标准迭代器循环。每个 1-D 循环完成后检查硬件错误标志。缓冲循环Buffered loop处理输入/输出数组未对齐或数据类型错误含字节交换而不符合底层 1-D 循环期望的情形数组假定为非连续的。其工作方式很像 strided-loop区别在于内层 1-D 循环被修改为按bufsize块对输入做预处理、对输出做后处理——bufsize是用户可设置参数。底层 1-D 计算循环作用于需要时复制过来的数据。该情形的设置与循环代码要复杂得多因为它必须处理临时缓冲区的内存分配决定输入输出数据是否使用缓冲区未对齐和/或类型错误时为需要缓冲区的输入/输出复制并可能转换数据对Object数组特判以便在需要复制/转换时正确处理引用计数把内层 1-D 循环拆分为bufsize大小的块可能有余数。同样每个 1-D 循环结束时检查硬件错误标志。最终输出处理Final output manipulationufunc 允许其他类数组类无缝通过接口某类的输入会诱导输出成为同类。机制如下若任一输入不是 ndarray 且定义了__array_wrap__方法则拥有最大__array_priority__属性的类决定所有输出的类型传入的输出数组除外。该输入数组的__array_wrap__方法会以 ufunc 返回的 ndarray 作为输入被调用。支持两种调用风格第一种ndarray 作为第一个参数元组形式的 context 作为第二个参数。context 为(ufunc, arguments, output argument number)——这是首先尝试的调用若抛出TypeError则改为只用 ndarray 作为第一个参数调用。方法Methodsufunc 有三种需要类似通用 ufunc 计算的方法ufunc.reduce、ufunc.accumulate和ufunc.reduceat。每个方法都需要一次设置 一次循环。方法有四种循环风格no-elements、one-element、strided-loop、buffered-loop——与通用函数调用的基本循环风格相同额外两种特例出现在输入数组分别有 0 个和 1 个元素时。设置三个方法的设置函数都是construct_reduce创建一个归约循环对象并填充完成循环所需的参数。所有方法只适用于2 输入、1 输出的 ufunc因此底层 1-D 循环按签名[otype, otype, otype]选择其中otype是请求的归约数据类型从每线程全局存储中取出缓冲大小和错误处理设置对小数组且未对齐或类型错误时制作副本以便使用非缓冲代码段选择循环策略数组有 1 或 0 个元素时选简单循环数组未对齐且类型正确则选 strided-loop否则必须使用 buffered-loop建立循环参数构造返回数组。输出数组形状依方法是reduce、accumulate还是reduceat而不同。若已提供输出数组则检查其形状若输出数组不是 C 连续、对齐且类型正确则制作带NPY_ARRAY_WRITEBACKIFCOPY标志的临时副本——方法可工作在规整的输出数组上函数完成调用PyArray_ResolveWritebackIfCopy时结果会复制回真正的输出数组最后设置迭代器以沿正确的轴依方法提供的 axis 参数循环然后返回实际计算例程。Reduce所有 ufunc 方法使用相同的底层 1-D 计算循环通过调整输入/输出参数实现相应归约。reduce 的关键1-D 循环被调用时输出与第二个输入指向内存中同一位置、步长均为 0第一个输入指向输入数组步长为所选轴的相应 stride。由此执行的操作是$$o i[0], \quad o i[k] \langle op \rangle o \quad (k1\ldots N)$$其中 $N1$ 是输入 $i$ 的元素数$o$ 是输出$i[k]$ 是 $i$ 沿所选轴的第 $k$ 个元素。对大于 1 维的数组该基本操作会重复使归约沿所选轴对每个 1-D 子数组进行——由移除所选维度的迭代器处理循环。缓冲循环必须小心调用循环函数前先复制并转换数据因为底层循环期望对齐的、类型含字节序正确的数据缓冲循环必须在不超过用户指定bufsize的块上处理复制与转换。Accumulateaccumulate与reduce非常相似输出与第二个输入都指向输出。区别在于第二个输入指向比当前输出指针落后一个 stride 的内存。因此执行的操作是$$o[0] i[0], \quad o[k] i[k] \langle op \rangle o[k-1] \quad (k1\ldots N)$$输出与输入形状相同当所选轴形状为 $N1$ 时每个 1-D 循环作用于 $N$ 个元素。同样缓冲循环在调用底层 1-D 计算循环前复制并转换数据。Reduceatreduceat是reduce与accumulate的推广它对输入数组中由索引指定的范围实现reduce。额外传入的 indices 参数会在循环计算前被检查确保每个输入沿所选维度都不超出输入数组范围。循环实现使用与reduce非常相似的代码重复执行indices 输入中元素个数次。特别地传给底层 1-D 计算循环的第一个输入指针指向索引数组指示的正确位置输出指针与第二个输入指针指向内存中同一位置1-D 计算循环的大小固定为当前索引与下一个索引之差当前索引为最后一个时下一个索引假定为数组沿所选维度的长度。由此 1-D 循环就在指定索引上实现一次reduce。未对齐或循环数据类型与输入/输出不匹配时使用缓冲代码处理数据复制到临时缓冲区必要时在调用底层 1-D 函数前转换为正确类型。临时缓冲区按元素大小创建不超过用户可设置的 buffer-size 值。因此循环必须足够灵活多次调用底层 1-D 计算循环以不超过 buffer-size 的块完成总计算。延伸阅读路径internals.code-explanations.rst本文的原始权威来源internals.rst 与 underthehood.rstNumPy 内部实现的更多讨论数据类型参考doc/source/reference/arrays.dtypes.rst、标量参考doc/source/reference/arrays.scalars.rst、迭代器参考doc/source/reference/arrays.nditer.rst、ufunc 参考doc/source/reference/ufuncs.rst核心 C 源码numpy/_core/src/multiarray/iterators.c迭代器与PyArray_Broadcast、numpy/_core/src/multiarray/mapping.c索引与PyArray_MapIterNew、numpy/_core/src/umath/ufunc_object.cufunc 调度与PyUFunc_GenericFunctionInternalC 代码风格约定可参考 doc/C_STYLE_GUIDE.rst。理解这些底层机制的价值在于strided 内存模型解释了为何切片视图零拷贝、NPY_ARRAY_ALIGNED解释了跨平台的内存安全边界、零 stride 广播解释了np.broadcast_to等操作为何几乎零开销而 ufunc 的三种循环策略则直接决定了np.add等操作在不同内存布局下的性能表现——这正是从会调用 NumPy进阶到理解 NumPy的关键一步。【免费下载链接】numpyThe fundamental package for scientific computing with Python.项目地址: https://gitcode.com/gh_mirrors/nu/numpy创作声明：本文部分内容由AI辅助生成（AIGC），仅供参考