如何将线性索引转换为支持负步长的下标

Question

有没有一种算法可以将线性索引转换成支持负步幅的下标列表？

出身背景

环境（如MATLAB、Julia等）和库（如NumPy）为跨步阵列（也称为NDArray）提供支持。跨步数组由线性内存（例如，单个底层缓冲区）支持，这与嵌套数组形成对比，其中每个嵌套数组对应一个维度。例如，考虑以下2x2矩阵

[ 1 2
  3 4 ]

实现为数组的数组

A = [ [ 1, 2 ], [ 3, 4 ] ]

在哪里（使用从零开始的索引）

a01 = A[0][1] = 2

我们可以将相同的 2x2 矩阵表示为跨频数组，如下所示（假设以行为主）

A = [ 1, 2,
      3, 4 ]

在哪里

a01 = A[ 2*0 + 1*1 ] = 2

通常，对于带条纹的NxM矩阵，元素（i， j）可以通过以下方式访问

function get( i, j ) {
    return buffer[ si*i + sj*j ];
}

其中缓冲区是底层数据缓冲区，si和sj分别对应于i和j维度上的步幅。假设一个行主步幅数组，对于上面的2x2矩阵，si=2和sj=1（省略元素字节长度）。

通常，可以从阵列形状计算跨距，如下所示：

function shape2strides( shape, order ) {
    var out = new Array( shape.length );
    var s = 1;
    var i;
    if ( order === 'column-major' ) {
        for ( i = 0; i < shape.length; i++ ) {
            out[ i ] = shape[ i ];
            s *= shape[ i ];
        }
        return out;
    } else { // row-major
        for ( i = shape.length-1; i >= 0; i-- ) {
            out[ i ] = shape[ i ];
            s *= shape[ i ];
        }
    }
}

为了方便使用跨步数组，环境/库通常提供方便的函数，允许在线性索引和下标之间轻松转换。例如，在MATLAB中，将下标转换为线性索引

idx = sub2ind( size( A ), i, j )

类似地，在MATLAB中从线性索引转换为下标

s = ind2sub( size( A ), idx )

Julia 也有 sub2ind 和 ind2sub。在NumPy中，您可以使用ravel_multi_index和unravel_index。

除了数据局部性之外，跨步数组也很方便，因为它们允许通过控制跨步是负还是正来创建数组“视图”。当步幅为负时，我们不是从左到右迭代，而是沿着那个维度从右到左迭代。为了支持这种迭代行为，我们需要确定底层数据缓冲区中第一个索引元素的位置。按照惯例，我们将这个索引称为“偏移量”，可以按如下方式计算

function strides2offset( shape, strides ) {
    var offset = 0;
    var i;
    for ( i = 0; i < shape.length; i++ ) {
        if ( strides[ i ] < 0 ) {
            offset -= strides[i] * ( shape[i]-1 ); // increments the offset
        }
    }
    return offset;
}

一旦我们有了偏移量，我们需要修改我们的get（ i， j ）函数，如下所示

function get( i, j ) {
    return buffer[ offset + si*i + sj*j ];
}

对于具有步幅2,1的2x2矩阵A，偏移量为0，从而返回上面的原始get函数。当步幅为2，-1时，偏移量为1；对于-2,1，偏移量为2；对于-2，-1，偏移量为3。因此，我们可以生成以下矩阵视图（假设行-主）

Dims: 2x2

Strides: 2,1
Offset: 0

A = [ 1, 2,
      3, 4 ]

Strides: 2,-1
Offset: 1

A = [ 2, 1,
      4, 3 ]

Strides: -2,1
Offset: 2

A = [ 3, 4,
      1, 2 ]

Strides: -2,-1
Offset: 3

A = [ 4, 3,
      2, 1 ]

以上视图突出了跨步阵列的优势之一：O（1）操作。例如，要将矩阵从左向右翻转，我们只需要翻转第二维度的步幅符号（假设行主要）。要上下翻转，我们翻转第一个维度的步幅符号（假设行主调）。要从左到右、从上到下翻转，我们翻转两个跨步的符号。所有上述操作不涉及触摸底层数据缓冲器；我们只是简单地改变了元数据的跨步数组。

sub2ind

从下标转换为线性索引很简单，即使考虑到负步幅（即步幅数组视图）。例如，对于任意维度的步幅数组，

function sub2ind( ...subscripts ) {
    var sub;
    var idx;
    var s;
    var n;

    idx = offset;
    for ( n = 0; n < shape.length; n++ ) {
        sub = subscripts[ n ];
        s = strides[ n ];
        if ( s < 0 && offset === 0 ) { // assume want "view" index
            idx -= sub * s; // always increments `idx`
        } else { // assume want underlying data buffer index
            idx += sub * s; // may increment or decrement `idx`
        }
    }
    return idx;
}

这里，我们允许从视图的角度或底层数据缓冲区的角度返回线性索引。当“offset”为< code>0时，我们假设我们总是将一个线性索引返回到视图中(它可能不对应于底层数据缓冲区中的线性索引)。换句话说，对于2x2矩阵视图，< code>(0，0)= 1

您可以向自己证明，在将元素查找扩展到负步幅时，sub2ind 为上述任何偏移返回相同的索引。

有关示例实现，请参阅 Julia、NumPy 和 stdlib。

ind2sub

这里提出的问题是如何实现sub2ind的反向，并支持负跨步。

对于正步幅（因此，偏移量为 0），我们可以使用模算术来恢复下标。例如，考虑用于解析 NxMxL 跨频阵列的线性索引的等式。

idx = offset + si*i + sj*j + sk*k

其中，假设行主要，si=nj*nk，sj=nk，sk=1和ni，nj，nk，分别是尺寸N，M，L。替换值，

idx = 0 + (nj*nk)*i + nk*j + k

可以重新排列

idx = nk*(nj*i + j) + k

如果我们使用nk取两边的模，

idx % nk = k

知道k，让我们重新排列初始方程

(idx - k) = nk*(nj*i + j)
(idx - k)/nk = nj*i + j

如果我们使用nj取两边的模，

((idx - k)/nk) % nj = j

知道了< code>j之后，让我们重新安排初始方程来求解< code>i

(((idx - k)/nk) - j)/nj = i

上述算法可以推广到任意数量的维度，并且易于实现（另请参阅Julia和NumPy）。

function ind2sub( idx, order ) {
    var out = new Array( shape.length );
    var s;
    var i;
    if ( order === 'column-major' ) {
        for ( i = 0; i < shape.length; i++ ) {
            s = idx % shape[ i ];
            idx -= s;
            idx /= shape[ i ];
            out[ i ] = s;
        }
    } else { // row-major
        for ( i = shape.length-1; i >= 0; i-- ) {
            s = idx % shape[ i ];
            idx -= s;
            idx /= shape[ i ];
            out[ i ] = s;
        }
    }
    return out;
}

然而，上面使用模算术的算法不支持负步幅。如果我们使用上面相同的过程来求解下标i， j， k，我们将从方程开始

idx = offset + nk*(nj*i + j) + k

可以简化为

idx-offset = nk*(nj*i + j) + k

当然，这里的问题是idx-偏移可能是负数，并有效地转移了可能的i， j， k值的范围（i应该在半开区间[0， N）；区间[0， M）上的j；和区间[0， L）上的k）。

这就提出了一个问题:是否存在一种算法，可以将线性索引转换成支持负步幅的下标。或者换句话说，有没有一种算法，给定底层数据缓冲区的线性索引，可以返回相应的视图下标？

其他语言/库中的实现(如Julia和NumPy)似乎只支持< code>offset = 0的情况。我在寻找一些更通用的东西，也可以应用于跨阵列视图。

任何指向现有实现/算法的指针都将受到极大的赞赏。

Answer 1

我写了一段代码，我认为可以解决你的问题。您要求的函数不止一个，因为其他函数与您的函数略有不同，因此我将它们全部放在这里以避免潜在的不一致。以下代码不是用任何特定语言编写的。但是，您可以在其中找到C语法的一些元素。

function calcStrides(shape[]) {
    strides[]; # Result array
    currentStride = 1;
    for(i = shape.size; 0 < i;) {
        --i;
        if(0 < shape[i]) {
            strides[i] = currentStride;
            currentStride *= shape[i];
        } else {
            strides[i] = -currentStride;
            currentStride *= -shape[i];
        }
    }
    return strides;
}

function calcOffset(shape[], strides[]) {
    offset = 0;
    for(i = 0; i < shape.size; ++i) {
        if(shape[i] < 0) {
            offset += strides[i] * (shape[i] + 1);
        }
    }
    return offset;
}

function sub2ind(strides[], offset, subs[]) {
    ind = offset;
    for(i = 0; i < strides.size; ++i) {
        ind += strides[i] * subs[i];
    }
    return ind;
};

function ind2sub(shape[], strides[], ind) {
    subs[]; # Result array
    for(i = 0; i < shape.size; ++i) {
        if(0 < strides[i]) {
            subs[i] = ind / strides[i];
            ind -= subs[i] * strides[i];
        } else {
            absSub = ind / -strides[i];
            subs[i] = -shape[i] - 1 - absSub;
            ind -= absSub * -strides[i];
        }
    }
    return subs;
}

Answer 2

(edit - I可能正在处理更简单的nd-index到flat index的情况，而您关注的是相反的情况。现在探究这个任务已经太晚了——我将在早上再次讨论这个问题。)

如果偏移量是正确的，我认为将n-d索引转换为平坦索引的相同公式适用于负和正步幅：

比较（3,4）数组的索引与其双翻转：

In [32]: x = np.arange(12).reshape(3,4)
In [33]: y = x[::-1, ::-1]
In [34]: x.strides
Out[34]: (16, 4)
In [35]: y.strides
Out[35]: (-16, -4)

可以使用< code>__array_interface__找到数据缓冲区“start”:

In [36]: x.__array_interface__['data']
Out[36]: (166934688, False)
In [37]: y.__array_interface__['data']
Out[37]: (166934732, False)
In [38]: 732-688
Out[38]: 44

缓冲区中有48个字节，但< code>y的偏移量是44，即x[2，3]元素的开始(本例中为11)。

现在测试＜code＞x＜code＞元素的平面索引：

In [39]: x[1,2]
Out[39]: 6            # value
In [40]: x[1:2, 2:3].__array_interface__['data']    # a view
Out[40]: (166934712, False)
In [41]: 688+(1*16)+(2*4)    # offset + si*i + sj*j 
Out[41]: 712

现在对y执行相同操作：

In [42]: y[1:2, 2:3].__array_interface__['data']
Out[42]: (166934708, False)
In [43]: 732+(1*-16)+(2*-4)
Out[43]: 708