使用涂料矢量来访问多维数组的任意轴向切片?

我正在构建一组函数来处理多维数组数据结构 ,我希望能够定义数组的任意切片 ,这样我就可以实现两个任意矩阵(又名Tensorsnd数组 )的广义内积。

我读过的一篇APL论文(我老实说找不到哪篇 – 我读过这么多篇幅)定义了左边矩阵X上的矩阵乘积,其尺寸为A;B;C;D;E;F和右矩阵Y尺寸G;H;I;J;K其中F==G as

 Z <- X +.× Y Z[A;B;C;D;E;H;I;J;K] <- +/ X[A;B;C;D;E;*] × Y[*;H;I;J;K] 

其中+/和的和 ,并且×将逐个元素应用于两个相同长度的向量。

所以我需要左边的“行”切片和右边的“列”切片。 我当然可以使用转置,然后使用“行”切片来模拟“列”切片,但我宁愿更优雅地做。

维基百科关于切片的文章引出了关于涂料载体的存根,这似乎是我正在寻找的奇迹治疗方法,但是没有太多可以继续下去。

如何使用涂料矢量来实现任意切片?

(很久以后我注意到Stride的一个数组有一些细节。)

定义

通过dope向量或描述符引用每个数组,可以实现通用数组切片(无论是否构建到语言中) – 包含第一个数组元素的地址的记录,然后是每个索引的范围和相应的系数。索引公式。 该技术还允许立即数组转置,索引反转,子采样等。对于像C这样的语言,其中索引总是从零开始,具有d索引的数组的掺杂矢量具有至少1 + 2d参数。
http://en.wikipedia.org/wiki/Array_slicing#Details

这是一个密集的段落,但它实际上都在那里。 所以我们需要一个像这样的数据结构:

 struct { TYPE *data; //address of first array element int rank; //number of dimensions int *dims; //size of each dimension int *weight; //corresponding coefficient in the indexing formula }; 

其中TYPE是元素类型的任何东西,矩阵的字段 。 为了简单和具体,我们只使用int 。 出于我自己的目的,我已经设计了各种类型的编码到整数句柄中,因此int做了工作,YMMV。

 typedef struct arr { int rank; int *dims; int *weight; int *data; } *arr; 

所有指针成员都可以在与结构本身相同的已分配内存块中分配位置(我们将其称为标题)。 但是通过替换早期使用的偏移和struct-hackery,算法的实现可以独立于块内(或没有)内的实际内存布局。

自包含数组对象的基本内存布局是

 rank dims weight data dims[0] dims[1] ... dims[rank-1] weight[0] weight[1] ... weight[rank-1] data[0] data[1] ... data[ product(dims)-1 ] 

间接数组共享数据(整个数组或1个或多个行切片)将具有以下内存布局

 rank dims weight data dims[0] dims[1] ... dims[rank-1] weight[0] weight[1] ... weight[rank-1] //no data! it's somewhere else! 

并且包含沿另一轴的正交切片的间接数组将具有与基本间接数组相同的布局,但具有适当修改的调光和重量。

索引(i0 i1 … iN)的元素的访问公式是

 a->data[ i0*a->weight[0] + i1*a->weight[1] + ... + iN*a->weight[N] ] 

假设每个索引i [j]在[0 … dims [j])之间。

在通常布局的行主arrays的weight向量中,每个元素是所有较低维度的乘积。

 for (int i=0; i 

因此,对于3×4×5arrays,元数据将是

 { .rank=3, .dims=(int[]){3,4,5}, .weight=(int[]){4*5, 5, 1} } 

或者对于7×6×5×4×3×2arrays,元数据将是

 { .rank=6, .dims={7,6,5,4,3,2}, .weight={720, 120, 24, 6, 2, 1} } 

施工

因此,要创建其中一个,我们需要与上一个问题相同的辅助函数来计算维度列表中的大小。

 /* multiply together rank integers in dims array */ int productdims(int rank, int *dims){ int i,z=1; for(i=0; i 

要分配,只需要malloc足够的内存并设置指向适当位置的指针。

 /* create array given rank and int[] dims */ arr arraya(int rank, int dims[]){ int datasz; int i; int x; arr z; datasz=productdims(rank,dims); z=malloc(sizeof(struct arr) + (rank+rank+datasz)*sizeof(int)); z->rank = rank; z->dims = z + 1; z->weight = z->dims + rank; z->data = z->weight + rank; memmove(z->dims,dims,rank*sizeof(int)); for(x=1, i=rank-1; i>=0; i--){ z->weight[i] = x; x *= z->dims[i]; } return z; } 

使用前一个答案中的相同技巧,我们可以创建一个变量参数接口,以简化使用。

 /* load rank integers from va_list into int[] dims */ void loaddimsv(int rank, int dims[], va_list ap){ int i; for (i=0; i 

甚至通过使用ppnarg的强大function计算其他参数来自动生成论证。

 #define array(...) (array)(PP_NARG(__VA_ARGS__),__VA_ARGS__) /* create a new array with specified dimensions */ 

现在构建其中一个非常容易。

 arr a = array(2,3,4); // create a dynamic [2][3][4] array 

访问元素

通过对elema的函数调用来检索元素,该函数将每个索引乘以相应的权重,对它们求和,并索引data指针。 我们返回一个指向元素的指针,以便调用者可以读取或修改它。

 /* access element of a indexed by int[] */ int *elema(arr a, int *ind){ int idx = 0; int i; for (i=0; irank; i++){ idx += ind[i] * a->weight[i]; } return a->data + idx; } 

相同的varargs技巧可以使一个更容易的接口elem(a,i,j,k)

轴向切片

首先,我们需要一种方法来指定要提取的维度和要折叠的维度。 如果我们只需要从维度中选择单个索引或所有元素,那么slice函数可以将rank int作为参数,并将-1解释为选择整个维度或0 ..将i -1作为选择单个索引。

 /* take a computed slice of a following spec[] instructions if spec[i] >= 0 and spec[i] < a->rank, then spec[i] selects that index from dimension i. if spec[i] == -1, then spec[i] selects the entire dimension i. */ arr slicea(arr a, int spec[]){ int i,j; int rank; for (i=0,rank=0; irank; i++) rank+=spec[i]==-1; int dims[rank]; int weight[rank]; for (i=0,j=0; i=a->rank) break; dims[i] = a->dims[j]; weight[i] = a->weight[j]; } arr z = casta(a->data, rank, dims); memcpy(z->weight,weight,rank*sizeof(int)); for (j=0; jrank; j++){ if (spec[j]!=-1) z->data += spec[j] * a->weight[j]; } return z; } 

收集与参数数组中的-1s对应的所有维度和权重,并用于创建新的数组头。 所有参数> = 0乘以它们的相关权重并添加到data指针,使指针偏移到正确的元素。

就数组访问公式而言,我们将其视为多项式。

 offset = constant + sum_i=0,n( weight[i] * index[i] ) 

因此,对于我们选择单个元素(+所有较低维度)的任何维度,我们将现在常量索引的因子分解并将其添加到公式中的常量项(在我们的C表示中是data指针本身)。

辅助函数casta使用共享data创建新的数组头。 slicea当然会改变权重值,但是通过计算权重本身, casta成为更普遍可用的函数。 它甚至可以用于构造一个动态数组结构,该结构直接在常规的C风格多维数组上运行,从而进行转换

 /* create an array header to access existing data in multidimensional layout */ arr casta(int *data, int rank, int dims[]){ int i,x; arr z=malloc(sizeof(struct arr) + (rank+rank)*sizeof(int)); z->rank = rank; z->dims = z + 1; z->weight = z->dims + rank; z->data = data; memmove(z->dims,dims,rank*sizeof(int)); for(x=1, i=rank-1; i>=0; i--){ z->weight[i] = x; x *= z->dims[i]; } return z; } 

转置

涂料矢量也可用于实现转座。 维度(和索引)的顺序可以更改。

请记住,这不像其他人那样正常的“转置”。 我们根本不重新排列数据。 这更恰当地称为'dope-vector pseudo-transpose'。 我们只需更改访问公式中的常量,重新排列多项式的系数,而不是更改数据。 从实际意义上讲,这只是对交换性和加法相关性的应用。

因此,为了具体,假设数据从假设地址500开始顺序排列。

 500: 0 501: 1 502: 2 503: 3 504: 4 505: 5 506: 6 

如果a是秩2,变暗{1,7),权重(7,1),则索引之和乘以添加到初始指针(500)的相关权重,产生每个元素的适当地址

 a[0][0] == *(500+0*7+0*1) a[0][1] == *(500+0*7+1*1) a[0][2] == *(500+0*7+2*1) a[0][3] == *(500+0*7+3*1) a[0][4] == *(500+0*7+4*1) a[0][5] == *(500+0*7+5*1) a[0][6] == *(500+0*7+6*1) 

因此,dope-vector伪转置重新排列权重和维度以匹配索引的新排序,但总和保持不变。 初始指针保持不变。 数据不会移动。

 b[0][0] == *(500+0*1+0*7) b[1][0] == *(500+1*1+0*7) b[2][0] == *(500+2*1+0*7) b[3][0] == *(500+3*1+0*7) b[4][0] == *(500+4*1+0*7) b[5][0] == *(500+5*1+0*7) b[6][0] == *(500+6*1+0*7) 

内在产品(又称矩阵乘法)

因此,通过使用通用切片或转置+“行” - 更复杂(更容易),可以实现广义内积。

首先,我们需要两个辅助函数,用于对产生矢量结果的两个矢量应用二元运算,并使用二元运算减少矢量,从而产生标量结果。

与上一个问题一样,我们将传入运算符,因此可以将许多不同的运算符使用相同的函数。 对于这里的样式,我将运算符作为单个字符传递,因此已经存在从C运算符到函数运算符的间接映射。 这是一个x-macro表 。

 #define OPERATORS(_) \ /* f F id */ \ _('+',+,0) \ _('*',*,1) \ _('=',==,1) \ /**/ #define binop(X,F,Y) (binop)(X,*#F,Y) arr (binop)(arr x, char f, arr y); /* perform binary operation F upon corresponding elements of vectors X and Y */ 

表中的额外元素用于null向量参数的reduce函数。 在这种情况下,reduce应返回运算符的标识元素 ,0表示+ ,1表示*

 #define reduce(F,X) (reduce)(*#F,X) int (reduce)(char f, arr a); /* perform binary operation F upon adjacent elements of vector X, right to left, reducing vector to a single value */ 

因此, binop会对操作符进行循环和切换。

 /* perform binary operation F upon corresponding elements of vectors X and Y */ #define BINOP(f,F,id) case f: *elem(z,i) = *elem(x,i) F *elem(y,i); break; arr (binop)(arr x, char f, arr y){ arr z=copy(x); int n=x->dims[0]; int i; for (i=0; i 

如果有足够的元素,则reduce函数会执行向后循环,如果有的话,将初始值设置为最后一个元素,并将初始值预设给运算符的identity元素。

 /* perform binary operation F upon adjacent elements of vector X, right to left, reducing vector to a single value */ #define REDID(f,F,id) case f: x = id; break; #define REDOP(f,F,id) case f: x = *elem(a,i) F x; break; int (reduce)(char f, arr a){ int n = a->dims[0]; int x; int i; switch(f){ OPERATORS(REDID) } if (n) { x=*elem(a,n-1); for (i=n-2;i>=0;i--){ switch(f){ OPERATORS(REDOP) } } } return x; } #undef REDID #undef REDOP 

通过这些工具,内部产品可以更高层次的方式实施。

 /* perform a (2D) matrix multiplication upon rows of x and columns of y using operations F and G. Z = X FG Y Z[i,j] = F/ X[i,*] G Y'[j,*] more generally, perform an inner product on arguments of compatible dimension. Z = X[A;B;C;D;E;F] +.* Y[G;H;I;J;K] |(F = G) Z[A;B;C;D;E;H;I;J;K] = +/ X[A;B;C;D;E;*] * Y[*;H;I;J;K] */ arr (matmul)(arr x, char f, char g, arr y){ int i,j; arr xdims = cast(x->dims,1,x->rank); arr ydims = cast(y->dims,1,y->rank); xdims->dims[0]--; ydims->dims[0]--; ydims->data++; arr z=arraya(x->rank+y->rank-2,catv(xdims,ydims)->data); int datasz = productdims(z->rank,z->dims); int k=y->dims[0]; arr xs = NULL; arr ys = NULL; for (i=0; irank+y->rank]; vector_index(i,z->dims,z->rank,idx); int *xdex=idx; int *ydex=idx+x->rank-1; memmove(ydex+1,ydex,y->rank); xdex[x->rank-1] = -1; free(xs); free(ys); xs = slicea(x,xdex); ys = slicea(y,ydex); z->data[i] = (reduce)(f,(binop)(xs,g,ys)); } free(xs); free(ys); free(xdims); free(ydims); return z; } 

此函数还使用了cast函数,它为casta提供了varargs接口。

 /* create an array header to access existing data in multidimensional layout */ arr cast(int *data, int rank, ...){ va_list ap; int dims[rank]; va_start(ap,rank); loaddimsv(rank,dims,ap); va_end(ap); return casta(data, rank, dims); } 

它还使用vector_index将1D索引转换为索引的nD向量。

 /* compute vector index list for ravel index ind */ int *vector_index(int ind, int *dims, int n, int *vec){ int i,t=ind, *z=vec; for (i=0; i 

github文件 。 其他支持function也在github文件中。


这个Q / A对是一系列相关问题的一部分,这些问题在实现包含用C语言编写的APL语言的解释器时出现。其他: 如何使用动态分配的任意维数组? ,以及如何将C数学运算符(+ - * /%)传递给函数result = mathfunc(x,+,y);? 。 其中一些材料以前发布到comp.lang.c。 comp.lang.apl中的更多背景信息。