NumPy 中 的 多 维 迭 代 器 代 码 之 美 第 19 章 江 涌
NumPy 是 Python 的 一 个 可 选 安 装 包, 提 供 了 一 个 功 能 强 大 的 N 维 数 组 对 象 NumPy 提 供 了 多 种 数 组 的 数 学 化 操 作 与 结 构 化 操 作, 使 得 Python 能 够 很 好 地 开 发 一 些 关 键 的 并 且 要 求 运 行 速 度 很 快 的 工 程 代 码 和 科 学 代 码 NumPy 中 通 过 切 片 (slicing) 的 概 念 来 实 现 快 速 结 构 化 操 作 记 法 为 [start:stop:stride], 例 如 im2=im[8:2:-1,9:1:-3] 按 照 slicing 方 式 选 取 的 新 影 像 将 与 原 始 影 像 共 享 数 据, 不 会 生 成 一 个 副 本, 减 少 计 算 机 资 源 的 消 耗
关 键 挑 战 经 常 需 要 遍 历 数 组 中 的 元 素, 在 遍 历 中 进 行 所 需 要 的 操 作 简 单 想 法 : 用 单 层 for 循 环 处 理 一 维, 双 层 for 循 环 处 理 两 维 但 当 维 数 N 是 一 个 任 意 整 数, 怎 么 办? 递 归 : 递 归 条 件 (recursive case), 基 线 条 件 (base case) Copy_ND(a,b,N) // 将 N 维 数 组 b 复 制 到 N 维 数 组 a 递 归 实 现 :if (N==0) copy memory from b to a return for i=0 to size of first dimension of a and b ptr_b=b[i] Copy_ND(ptr_a,ptr_b,N-1) a[i]=ptr_a
递 归 算 法 在 每 次 迭 代 中 进 行 函 数 调 用, 容 易 产 生 速 度 很 慢 的 代 码 ; 许 多 算 法 需 要 保 存 中 间 值 用 于 后 续 的 递 归 调 用 ( 求 最 大 值 ), 这 些 值 将 被 作 为 递 归 调 用 的 参 数 传 递, 很 难 提 供 用 于 递 归 解 决 方 案 的 简 化 工 具 因 此,NumPy 使 用 迭 代 来 完 成 迭 代 器 (Iterator) 是 一 种 简 化 这 些 算 法 的 抽 象, 包 含 了 单 个 循 环 内 遍 历 数 组 中 所 有 元 素 的 思 想 迭 代 器 的 两 个 基 本 方 法 : hasnext 是 否 还 有 下 一 个 元 素 ; next 返 回 下 一 个 元 素 for x in iterobj: process(x)
数 组 的 内 存 模 型 邻 接 型 数 组 : 在 内 存 中 连 续 存 放 一 个 二 维 4*5 数 组 >>>p=[[1,2,3,4,5], [6,7,8,9,10], [11,12,13,14,15], [16,17,18,19,20]] >>>from numpy import * >>>pp=array(p) >>>p1=pp[1:3,1:4] 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
跨 度 数 组 中 某 一 维 的 跨 度 (stride): 沿 着 这 一 维, 或 者 说 数 轴, 从 数 组 中 的 一 个 元 素 移 动 到 下 一 个 元 素, 需 要 跳 过 多 少 字 节 ( 跨 度 可 以 是 负 数 ) pp: 第 一 维 跨 度 4*5, 第 二 维 跨 度 4; p1: 第 一 维 跨 度 4*5, 第 二 维 跨 度 4.
迭 代 器 设 计 迭 代 器 循 环 伪 码 : set up iterator (including pointing the current value to the first value in the array) while iterator not done: process the current value point the current value to the next value 设 计 分 为 三 部 分 : 1. Moving to the next value 2. Termination 3. Setup
递 进 确 定 按 怎 样 的 顺 序 来 提 取 元 素 NumPy 中 通 过 使 用 一 组 数 字 来 模 拟 简 单 计 数 而 实 现 用 N 元 整 数 组 来 表 示 当 前 位 置, 数 组 的 形 状 n1 n2... nn (0,,0) 表 示 数 组 第 一 个 元 素,( n1 1, n2 1,..., n N 1) 表 示 数 组 的 最 后 一 个 元 素 下 一 个 元 素 的 位 置 是 将 最 后 一 个 数 字 加 1 来 得 到 若 第 i 个 数 字 到 达 了 n i, 那 么 这 个 数 字 将 被 设 置 为 0, 而 第 (i-1) 个 数 字 将 增 1 比 如 对 于 数 组 来 说 3 2 3 (0,0,0)(0,0,1)(0,0,2)(0,1,0)(0,1,1)(0,1,2)(1,0,0)...(2,1,1)(2,1,2)
假 设 data 是 指 向 数 组 起 始 位 置 的 指 针,counter[N] 是 计 数 器 数 组,strides[N] 是 跨 度 值 数 组, 那 么 下 面 的 运 算 将 把 dataptr 设 置 为 指 向 数 组 当 前 值 的 第 一 个 字 节 : dataptr = (char *)data; for (i=0; i<n; i++) dataptr += counter[i]*strides[i]; 事 实 上 可 以 在 记 录 计 数 器 的 同 时 也 记 录 指 针, 当 计 数 器 的 第 i 个 下 标 增 1 时,dataptr 将 会 增 加 strides[] i ; 当 第 i 个 下 标 复 位 为 0 时, 数 组 当 前 值 的 内 存 地 址 应 该 减 去 ( n 1) strides[ i] i 迭 代 器 负 责 维 护 计 数 器 与 指 向 当 前 值 的 指 针
终 止 判 断 迭 代 器 何 时 完 成 及 如 何 发 出 终 止 信 号 1. 附 加 标 志 变 量 到 迭 代 器 上, 每 次 迭 代 中 都 进 行 判 断, 如 果 没 有 元 素 了 就 设 置 这 个 标 志 2. 查 找 在 第 一 维 的 计 数 器 中 从 ni 1 到 0 的 跃 迁 点 3. 如 果 给 定 了 数 组 的 大 小, 只 需 记 住 将 要 进 行 的 迭 代 次 数 NumPy 中 将 迭 代 的 总 次 数 作 为 信 息 保 存 下 来, 以 及 保 存 一 个 到 目 前 为 止 的 迭 代 次 数 的 动 态 计 数 器, 当 达 到 迭 代 的 总 次 数 时 迭 代 终 止
构 建 需 要 保 存 的 信 息 : 整 数 计 数 器 ( 初 始 设 置 为 0) 下 标 计 数 器 ( 初 始 设 置 为 (0,0,,0)) 判 断 是 否 基 于 简 单 的 邻 接 内 存, 设 置 标 志 保 存 判 断 结 果 为 了 加 速 回 卷 步 骤, 保 存 每 一 维 ( ni 1), 避 免 重 复 计 算 ( ni 1) strides[ i] 跨 度 信 息, 维 数 信 息, 元 素 数 量 保 存 系 数 li ( i = 1,... N) 来 简 化 整 数 计 数 器 与 下 标 计 数 器 的 转 化 当 前 指 针
数 组 中 的 每 一 项 都 可 以 用 一 个 在 0 和 n1... n N 之 间 1 的 整 数 k 或 者 下 标 计 数 器 ( k1,..., k N ) 来 表 示, 这 个 关 系 可 以 被 定 义 为 N N ( ) k k n = i i= 1 j=+ i 1 j 或 l1 = k N l mod( ) i = l i 1 nj k i j= i i = N n j j = i + 1 l
coords[n] dims_m1[n] strides[n] backstrides[n] nd_m1 dataptr size index factors 下 标 计 数 器 N 维 数 组 N 维 数 组 保 存 每 一 维 ni 1 N 维 数 组, 保 存 每 维 跨 度 N 维 数 组, 回 卷 时 需 要 移 动 的 数 量 strides[i] 维 数 ( ni 1) 指 向 当 前 位 置 内 存 的 指 针 所 有 元 素 数 量 n n n 1 2... N 整 数 计 数 器 从 0 变 到 size-1 计 算 一 维 下 标 和 N 维 下 标 转 换 时 辅 助 数 组
记 录 迭 代 计 数 器 注 意 计 数 器 的 递 进 都 是 从 最 后 一 维 增 1 开 始, 当 某 维 大 于 n 1 时 发 生 回 卷, 此 时 可 能 使 其 他 维 的 下 标 也 发 生 回 卷 i 于 是 可 以 从 最 后 维 开 始 向 前 循 环 在 当 前 维 上 判 断 当 前 下 标 是 不 是 小 于 n 1, 如 果 是, 则 将 下 标 位 置 加 1, 并 且 将 当 前 维 对 应 的 strides[i] 加 到 dataptr 上, 跳 出 循 环 ; 若 第 i 维 下 标 的 计 数 器 大 于 或 等 于 i n 1 i, 就 重 新 设 为 0, 且 将 dataptr 减 去 该 维 对 应 的 backstrides[i]( 回 卷 ), 继 续 循 环 判 断 前 一 维
C 语 言 实 现 for (i=it->nd_m1; i>=0; i--) { if (it->coords[i] < it->dims_m1[i]) { it->coords[i]++; it->dataptr += it->strides[i]; break; } else { it->coords[i] = 0; it->dataptr -= it->backstrides[i]; } }
使 用 while 语 句 done = 0; i = it->nd_ml; while (!done i>=0) { /*&&*/ if (it->coords[i] < it->dims_m1[i]) { it->coords[i]++; it->dataptr += it->strides[i]; done = 1; } else { it->coords[i] = 0; it->dataptr -= it->backstrides[i];} i--; }
NumPy 迭 代 器 的 结 构 typedef struct { PyObject_HEAD int nd_m1; npy_intp index, size; npy_intp coords[npy_maxdims]; npy_intp dims_m1[npy_maxdims]; npy_intp strides[npy_maxdims]; npy_intp backstrides[npy_maxdims]; npy_intp factors[npy_maxdims]; PyArrayObject *ao; char *dataptr; npy_bool contiguous; } PyArrayIterObject; 指 向 构 建 迭 代 器 的 原 始 数 组 的 指 针 判 断 是 否 邻 接 数 组
接 口 在 NumPy 中 it=pyarray_iternew(ao) ------------- 构 建 数 组 ao 的 迭 代 器 PyArray_ITER_NOTDONE(it) ------------- 判 断 迭 代 是 否 结 束 PyArray_ITER_NEXT(it) ------------- 实 现 迭 代 的 下 一 个 位 置 PyArray_ITER_DATA(it) ------------- 得 到 指 向 当 前 值 第 一 个 字 节 的 指 针
示 例 : 计 算 N 维 数 组 中 的 最 大 值 ( 假 设 数 组 ao 是 double 类 型 ) #include <float.h> double *currval, maxval=-dbl_max; PyArrayIterObject *it; it = PyArray_IterNew(ao); while (PyArray_ITER_NOTDONE(it)) { currval = (double *)PyArray_ITER_DATA(it); if (*currval > maxval) maxval = *currval; PyArray_ITER_NEXT(it); }
用 迭 代 器 处 理 邻 接 数 组 也 是 很 快 的, 但 是 更 快 的 还 是 传 统 办 法 : double *currval, maxval=-max_double; int size; currval = (double *)PyArray_DATA(ao); size = PyArray_SIZE(ao); while (size--) { if (*currval > maxval) maxval = *currval; currval += 1; }
迭 代 器 的 使 用 排 除 某 一 维 的 迭 代 当 操 作 不 涉 及 到 某 维 时 可 以 使 迭 代 器 跨 过 这 维 进 行 迭 代, 以 获 得 速 度 提 升 例 如 数 组 a[3][3], 只 想 修 改 a(0,0),a(1,0),a(2,0) 的 值 1 2 3 1 2 3
通 常 是 排 除 最 后 一 维,NumPy 的 前 身 Numeric 就 是 引 入 这 种 方 法 实 现 数 学 功 能 NumPy 中 对 迭 代 器 稍 微 改 动 it=pyarray_iterallbutaxia(array,&dim) dim 为 排 除 的 维 当 输 入 的 维 是 -1 时, 将 会 自 行 选 择 最 小 非 零 跨 度 值 的 维 另 一 种 经 常 的 选 择 就 是 排 除 有 最 大 元 素 数 量 的 维, 这 样 将 迭 代 次 数 降 至 最 低
实 现 : 1. 将 迭 代 大 小 除 以 将 要 移 除 的 维 的 长 度 ; 2. 将 被 选 择 的 维 的 元 素 数 量 设 置 为 1:dims_m1[i]=0; 3. 将 该 维 在 backstrides 中 相 应 位 置 上 的 值 设 置 为 0; 4. 将 邻 接 标 志 重 新 设 置 为 0, 因 为 现 在 要 处 理 的 数 组 在 内 存 中 将 不 会 是 邻 接 的 在 循 环 的 每 次 迭 代 中, 迭 代 器 将 指 向 数 组 所 选 择 维 的 第 一 个 元 素
迭 代 器 的 使 用 多 重 迭 代 数 组 加 法 需 要 多 个 迭 代 器, 每 个 输 入 数 组 一 个 迭 代 器, 并 且 输 出 数 组 一 个 迭 代 器 NumPy 提 供 了 一 个 多 重 迭 代 器, 可 以 对 多 个 迭 代 器 同 时 处 理 这 种 多 迭 代 器 经 常 用 于 自 动 处 理 NumPy 的 广 播 (broadcasting) 功 能 广 播 能 使 一 个 (4,1) 形 状 的 数 组 加 到 一 个 (3) 形 状 的 数 组 上, 得 到 形 状 (4,3) 的 数 组, 也 能 对 一 个 (5,1,1) 形 状 数 组, 一 个 (4,1) 形 状 数 组 和 一 个 (3) 形 状 数 组 进 行 运 算, 得 到 一 个 (5,4,3) 形 状 的 数 组
>>>p=array([[1,2,3], >>>a=array([[1],[2],[3],[4]]) [4,5,6], >>>a1=[2,1,4] [7,8,9]]) >>>print a+a1 >>>p1=[2,1,4] [[3,2,5], >>>print p+p1 [4,3,6], [[3,3,7], [5,4,7], [6,6,10], [6,5,8]] [9,9,13]]
The rules of broadcasting are: Arrays with fewer dimensions are treated as occupying the last dimensions of an array that has the full number of dimensions, so that all arrays have the same number of dimensions. The new, initial dimensions are filled in with 1s. The length of each dimension in the final broadcast shape is the greatest length of that dimension in any of the arrays.
For each dimension, all inputs must either have the same number of elements as the broadcast result or a 1 as the number of elements. Arrays with a single element in a particular dimension act as if that element were virtually copied to all positions during the iteration. In effect, the element is broadcast to the additional positions.
修 改 迭 代 器 实 现 广 播 功 能 : 修 改 迭 代 器 的 形 状 为 匹 配 广 播 的 形 状, 用 于 广 播 维 数 的 strides 和 backstrides 被 修 改 为 0, 这 样 这 一 维 递 进 时, 迭 代 器 不 会 移 动 数 据 指 针, 无 需 复 制 内 存
PyObject *multi; PyObject *in1, *in2; double *i1p, *i2p, *op; /* get in1 and in2 (assumed to be arrays of NPY_DOUBLE) */ /* first argument is the number of input arrays; the next (variable number of) arguments are the array objects */ multi = PyArray_MultiNew(2, in1, in2); /* construct output array */ out = PyArray_SimpleNew(PyArray_MultiIter_NDIM(multi), 输 出 数 组 PyArray_MultiIter_DIMS(multi), NPY_DOUBLE); 数 组 当 前 数 据 指 向 op = PyArray_DATA(out); 多 重 迭 代 器
while(pyarray_multiiter_notdone(multi)) { /* get (pointers to) the current value in each array */ i1p = PyArray_MultiIter_DATA(multi, 0); i2p = PyArray_MultiIter_DATA(multi, 1); /* perform the operation for this element */ *op = *ip1 + *ip2 op += 1; /* Advance output array pointer */ /* Advance all the input iterators */ PyArray_MultiIter_NEXT(multi); }
总 结 迭 代 器 是 一 个 非 常 漂 亮 的 抽 象 NumPy 的 迭 代 器 为 Python 的 数 组 运 算 提 供 了 很 强 大 的 灵 活 工 具, 并 且 无 需 考 虑 数 组 在 内 存 中 是 否 邻 接, 以 切 片 的 概 念 减 少 内 存 的 消 耗 它 的 优 化 循 环 与 广 播 机 制 都 是 其 漂 亮 的 发 光 点