使用掺杂向量访问多维数组的任意轴向切片？

2024-04-09

我正在构建一套功能来与多维数组数据结构 https://stackoverflow.com/questions/30023867/how-can-i-work-with-dynamically-allocated-arbitrary-dimensional-arrays/30023868#30023868我希望能够定义任意slices的数组，这样我就可以实现两个任意矩阵的广义内积（又名Tensors or n 维数组).

我读过的一篇 APL 论文（老实说我找不到哪篇——我读了很多篇）定义了左矩阵上的矩阵乘积X有尺寸A;B;C;D;E;F和右矩阵Y有尺寸G;H;I;J;K where 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]

where +/ is sum of、× 将逐个元素应用于两个长度相同的向量。

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

维基百科的文章关于slicing http://en.wikipedia.org/wiki/Array_slicing导致一个关于涂料矢量图 http://en.wikipedia.org/wiki/Dope_vector这似乎是我正在寻找的奇迹疗法，但没有太多可以继续的地方。

如何使用涂料向量实现任意切片？

（很久以后我才注意到数组的步幅 http://en.wikipedia.org/wiki/Stride_of_an_array其中有一些细节。）

定义

一般的数组切片可以通过通过 dope 向量或描述符引用每个数组来实现（无论是否内置在语言中）——一条记录包含第一个数组元素的地址，然后是每个索引的范围和相应的系数索引公式。该技术还允许立即数组转置、索引反转、二次采样等。对于像 C 这样的语言，索引总是从零开始，具有 d 个索引的数组的掺杂向量至少有 1 + 2d 个参数。
http://en.wikipedia.org/wiki/Array_slicing#Details 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
};

Where TYPE无论元素类型是什么，field矩阵的。为了简单和具体，我们将使用int。为了我自己的目的，我设计了一种各种类型的编码整数句柄 https://groups.google.com/d/topic/comp.lang.apl/KSwo40-eLKI/discussion so int做这份工作me, 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] ) 之间。

In the weight正常布局的向量行优先数组，每个元素都是所有较低维度的乘积。

for (int i=0; i<rank; i++)
    weight[i] = product(dims[i+1 .. rank-1]);

因此对于 3×4×5 数组，元数据将是

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

或者对于 7×6×5×4×3×2 数组，元数据将为

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

建造

因此，要创建其中之一，我们需要来自上一个问题 https://stackoverflow.com/questions/30023867/how-can-i-work-with-dynamically-allocated-arbitrary-dimensional-arrays从维度列表计算大小。

/* multiply together rank integers in dims array */
int productdims(int rank, int *dims){
    int i,z=1;
    for(i=0; i<rank; i++)
        z *= dims[i];
    return z;
}

分配的话，简单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<rank; i++){
        dims[i]=va_arg(ap,int);
    }
}

/* create a new array with specified rank and dimensions */
arr (array)(int rank, ...){
    va_list ap;
    //int *dims=calloc(rank,sizeof(int));
    int dims[rank];
    int i;
    int x;
    arr z;

    va_start(ap,rank);
    loaddimsv(rank,dims,ap);
    va_end(ap);

    z = arraya(rank,dims);
    //free(dims);
    return z;
}

甚至自动生成rank通过使用以下强大的力量来计算其他参数来计算参数ppnarg https://stackoverflow.com/a/6722112/733077.

#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; i<a->rank; i++){
        idx += ind[i] * a->weight[i];
    }
    return a->data + idx;
}

同样的可变参数技巧可以使界面更简单elem(a,i,j,k).

轴向切片

To take a slice, first we need a way of specifying which dimensions to extract and which to collapse. If we just need to select a single index or all elements from a dimension, then the slice function can take rank ints as arguments and interpret -1 as selecting the whole dimension or 0..dims_i-1 as selecting a single index.

/* 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; i<a->rank; i++)
        rank+=spec[i]==-1;
    int dims[rank];
    int weight[rank];
    for (i=0,j=0; i<rank; i++,j++){
        while (spec[j]!=-1) j++;
        if (j>=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; j<a->rank; j++){
        if (spec[j]!=-1)
            z->data += spec[j] * a->weight[j];
    }   
    return z;
}

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

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

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

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

辅助函数casta使用共享创建新的数组头data. slicea当然会改变权重值，但是通过计算权重本身，casta成为更普遍可用的功能。它甚至可以用来构造一个动态数组结构，直接对常规 C 风格的多维数组进行操作，从而casting.

/* 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;
}

转置

掺杂向量也可用于实现转置。维度（和索引）的顺序可以更改。

请记住，这不是像其他人一样的正常“转置” 做。我们根本不重新排列数据。这更正确地称为“掺杂向量伪转置”。我们不改变数据，只改变访问公式中的常数，重新排列多项式的系数。从真正意义上来说，这只是交换律的一个应用加法的结合性。

因此，为了具体起见，假设数据已排列按顺序从假设地址 500 开始。

如果 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)

所以涂料向量伪转置重新排列了重量和尺寸以匹配新的订购指数，但总和保持不变。最初的指针保持不变。数据不动。

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)

内积（又名矩阵乘法）

因此，通过使用通用切片或转置+“行”-切片（更容易），可以实现广义内积。

首先，我们需要两个辅助函数，用于将二元运算应用于两个向量，生成向量结果，并通过二元运算减少向量，生成标量结果。

正如在上一个问题 https://stackoverflow.com/questions/29960273/how-to-pass-a-c-math-operator-into-a-function-result-mathfuncx-y我们将传入运算符，因此同一函数可以与许多不同的运算符一起使用。对于此处的样式，我将运算符作为单个字符传递，因此已经存在从 C 运算符到我们的函数的运营商。这是一x宏表 https://stackoverflow.com/questions/6635851/real-world-use-of-x-macros/6636596#6636596.

#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 */

表中的额外元素用于reduce用于空向量参数情况的函数。在这种情况下，reduce 应该返回算子的单位元, 0 for +, 1 for *.

#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 */

So the 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<n; i++){
        switch(f){
            OPERATORS(BINOP)
        }
    }
    return z;
}
#undef BINOP

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

/* 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 F.G 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; i<datasz; i++){
        int idx[x->rank+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它提供了一个 varargs 接口casta.

/* 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<n; i++){
        z[n-1-i] = t % dims[n-1-i];
        t /= dims[n-1-i];
    }
    return z;
}

github文件 https://github.com/luser-dr00g/inca/blob/master/arrind.c#L404。其他支持功能也在 github 文件中。

该问答对是在实施过程中出现的一系列相关问题的一部分inca https://github.com/luser-dr00g/inca/用 C 编写的 APL 语言解释器。其他：如何使用动态分配的任意维数组？ https://stackoverflow.com/questions/30023867/how-can-i-work-with-dynamically-allocated-arbitrary-dimensional-arrays , and 如何将 C 数学运算符 (+-*/%) 传递到函数 result=mathfunc(x,+,y); 中？ https://stackoverflow.com/questions/29960273/how-to-pass-a-c-math-operator-into-a-function-result-mathfuncx-y。其中一些材料之前已发布到编译语言 https://groups.google.com/d/topic/comp.lang.c/Z0mycsYvPyI/discussion。更多背景编译语言.apl https://groups.google.com/d/msg/comp.lang.apl/zCPJ8KaZyWg/0U_pC9skLbUJ.

本文内容由网友自发贡献，版权归原作者所有，本站不承担相应法律责任。如您发现有涉嫌抄袭侵权的内容，请联系:hwhale#tublm.com(使用前将#替换为@)