Microsoft PowerPoint - seminar-10-ch19-iterator.ppt

Size: px

Start display at page:

Download "Microsoft PowerPoint - seminar-10-ch19-iterator.ppt"

纸孔
7 years ago
Views:

1 NumPy 中的多维迭代器代码之美第 19 章江涌

2 NumPy 是 Python 的一个可选安装包, 提供了一个功能强大的 N 维数组对象 NumPy 提供了多种数组的数学化操作与结构化操作, 使得 Python 能够很好地开发一些关键的并且要求运行速度很快的工程代码和科学代码 NumPy 中通过切片 (slicing) 的概念来实现快速结构化操作记法为 [start:stop:stride], 例如 im2=im[8:2:-1,9:1:-3] 按照 slicing 方式选取的新影像将与原始影像共享数据, 不会生成一个副本, 减少计算机资源的消耗

NumPy 中通过切片 (slicing) 的概念来实现快速结构化操作记法为 [start:stop:stride], 例如

3 关键挑战经常需要遍历数组中的元素, 在遍历中进行所需要的操作简单想法 : 用单层 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

递归 : 递归条件 (recursive case), 基线条件 (base case) Copy_ND(a,b,N) // 将 N 维数组 b 复制到 N 维数组

4 递归算法在每次迭代中进行函数调用, 容易产生速度很慢的代码 ; 许多算法需要保存中间值用于后续的递归调用 ( 求最大值 ), 这些值将被作为递归调用的参数传递, 很难提供用于递归解决方案的简化工具因此,NumPy 使用迭代来完成迭代器 (Iterator) 是一种简化这些算法的抽象, 包含了单个循环内遍历数组中所有元素的思想迭代器的两个基本方法 : hasnext 是否还有下一个元素 ; next 返回下一个元素 for x in iterobj: process(x)

用迭代来完成迭代器 (Iterator) 是一种简化这些算法的抽象, 包含了单个循环内遍历数组中所有元素的思想

5 数组的内存模型邻接型数组 : 在内存中连续存放一个二维 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]

>>>from numpy import * >>>pp=array(p) >>>p1=pp[1:3,1:4] 1 2 3 4 5 6 7

6 跨度数组中某一维的跨度 (stride): 沿着这一维, 或者说数轴, 从数组中的一个元素移动到下一个元素, 需要跳过多少字节 ( 跨度可以是负数 ) pp: 第一维跨度 4*5, 第二维跨度 4; p1: 第一维跨度 4*5, 第二维跨度 4.

7 迭代器设计迭代器循环伪码 : 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

done: process the current value point the current value to the next

8 递进确定按怎样的顺序来提取元素 NumPy 中通过使用一组数字来模拟简单计数而实现用 N 元整数组来表示当前位置, 数组的形状 n1 n2... nn (0,,0) 表示数组第一个元素,( n1 1, n2 1,..., n N 1) 表示数组的最后一个元素下一个元素的位置是将最后一个数字加 1 来得到若第 i 个数字到达了 n i, 那么这个数字将被设置为 0, 而第 (i-1) 个数字将增 1 比如对于数组来说 (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)

.., n N 1) 表示数组的最后一个元素下一个元素的位置是将最后一个数字加 1 来得到若第 i 个数字到达了 n i,

9 假设 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 迭代器负责维护计数器与指向当前值的指针

counter[i]*strides[i]; 事实上可以在记录计数器的同时也记录指针, 当计数器的第 i 个下标增 1 时,dataptr 将会增加

10 终止判断迭代器何时完成及如何发出终止信号 1. 附加标志变量到迭代器上, 每次迭代中都进行判断, 如果没有元素了就设置这个标志 2. 查找在第一维的计数器中从 ni 1 到 0 的跃迁点 3. 如果给定了数组的大小, 只需记住将要进行的迭代次数 NumPy 中将迭代的总次数作为信息保存下来, 以及保存一个到目前为止的迭代次数的动态计数器, 当达到迭代的总次数时迭代终止

11 构建需要保存的信息 : 整数计数器 ( 初始设置为 0) 下标计数器 ( 初始设置为 (0,0,,0)) 判断是否基于简单的邻接内存, 设置标志保存判断结果为了加速回卷步骤, 保存每一维 ( ni 1), 避免重复计算 ( ni 1) strides[ i] 跨度信息, 维数信息, 元素数量保存系数 li ( i = 1,... N) 来简化整数计数器与下标计数器的转化当前指针

12 数组中的每一项都可以用一个在 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

13 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 N 整数计数器从 0 变到 size-1 计算一维下标和 N 维下标转换时辅助数组

组, 回卷时需要移动的数量 strides[i] 维数 ( ni 1) 指向当前位置内存的指针所有元素

14 记录迭代计数器注意计数器的递进都是从最后一维增 1 开始, 当某维大于 n 1 时发生回卷, 此时可能使其他维的下标也发生回卷 i 于是可以从最后维开始向前循环在当前维上判断当前下标是不是小于 n 1, 如果是, 则将下标位置加 1, 并且将当前维对应的 strides[i] 加到 dataptr 上, 跳出循环 ; 若第 i 维下标的计数器大于或等于 i n 1 i, 就重新设为 0, 且将 dataptr 减去该维对应的 backstrides[i]( 回卷 ), 继续循环判断前一维

下标位置加 1, 并且将当前维对应的 strides[i] 加到 dataptr 上, 跳出循环 ; 若第 i 维下标的计数器大于或

15 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]; } }

16 使用 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--; }

17 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; 指向构建迭代器的原始数组的指针判断是否邻接数组

npy_intp backstrides[npy_maxdims]; npy_intp factors[npy_maxdims]; PyArrayObject *ao; char

18 接口在 NumPy 中 it=pyarray_iternew(ao) 构建数组 ao 的迭代器 PyArray_ITER_NOTDONE(it) 判断迭代是否结束 PyArray_ITER_NEXT(it) 实现迭代的下一个位置 PyArray_ITER_DATA(it) 得到指向当前值第一个字节的指针

19 示例 : 计算 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); }

$PyArray_IterNew(ao); while (PyArray_ITER_NOTDONE(it)) { currval = (double$

20 用迭代器处理邻接数组也是很快的, 但是更快的还是传统办法 : 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; }

21 迭代器的使用排除某一维的迭代当操作不涉及到某维时可以使迭代器跨过这维进行迭代, 以获得速度提升例如数组 a[3][3], 只想修改 a(0,0),a(1,0),a(2,0) 的值

22 通常是排除最后一维,NumPy 的前身 Numeric 就是引入这种方法实现数学功能 NumPy 中对迭代器稍微改动 it=pyarray_iterallbutaxia(array,&dim) dim 为排除的维当输入的维是 -1 时, 将会自行选择最小非零跨度值的维另一种经常的选择就是排除有最大元素数量的维, 这样将迭代次数降至最低

23 实现 : 1. 将迭代大小除以将要移除的维的长度 ; 2. 将被选择的维的元素数量设置为 1:dims_m1[i]=0; 3. 将该维在 backstrides 中相应位置上的值设置为 0; 4. 将邻接标志重新设置为 0, 因为现在要处理的数组在内存中将不会是邻接的在循环的每次迭代中, 迭代器将指向数组所选择维的第一个元素

24 迭代器的使用多重迭代数组加法需要多个迭代器, 每个输入数组一个迭代器, 并且输出数组一个迭代器 NumPy 提供了一个多重迭代器, 可以对多个迭代器同时处理这种多迭代器经常用于自动处理 NumPy 的广播 (broadcasting) 功能广播能使一个 (4,1) 形状的数组加到一个 (3) 形状的数组上, 得到形状 (4,3) 的数组, 也能对一个 (5,1,1) 形状数组, 一个 (4,1) 形状数组和一个 (3) 形状数组进行运算, 得到一个 (5,4,3) 形状的数组

25 >>>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]]

26 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.

27 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.

28 修改迭代器实现广播功能 : 修改迭代器的形状为匹配广播的形状, 用于广播维数的 strides 和 backstrides 被修改为 0, 这样这一维递进时, 迭代器不会移动数据指针, 无需复制内存

29 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); 多重迭代器

30 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); }

31 总结迭代器是一个非常漂亮的抽象 NumPy 的迭代器为 Python 的数组运算提供了很强大的灵活工具, 并且无需考虑数组在内存中是否邻接, 以切片的概念减少内存的消耗它的优化循环与广播机制都是其漂亮的发光点

《C语言基础入门》课程教学大纲

《C语言基础入门》课程教学大纲 C 语言开发入门教程课程教学大纲课程编号 :201409210011 学分 :5 学分学时 :58 学时 ( 其中 : 讲课学时 :39 学时上机学时 :19 学时 ) 先修课程 : 计算机导论后续课程 :C++ 程序设计适用专业 : 信息及其计算机相关专业开课部门 : 计