关于 c ：遍历任何 boost::multi_array 的除第 d 维之外的所有维度

Iterate over all but d-th dimension of any boost::multi_array

很多时候，人们想要沿着 N 维数组 A 的维 d 应用操作 f()。这意味着遍历 A 的所有剩余维度。我试图弄清楚 boost::multi_array 是否能够做到这一点。函数 f(A) 应该适用于所有类型的 boost::multi_array，包括 boost:multi_array_ref、boost::detail::multi_array::sub_array 和 boost::detail::multi_array::array_view，理想情况下也适用于右值类型，例如 boost::multi_array_ref<T, NDims>::reference。

我能想到的最好的方法是实现一个 reshape() 函数，该函数可用于将 ND 数组重塑为 3D 数组，这样工作维度始终是中间维度。这是 f.hpp:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57

#include”boost/multi_array.hpp”
#include <ostream>

using namespace boost;

typedef multi_array_types::index index_t;
typedef multi_array_types::index_range range;

template <template <typename, std::size_t, typename…> class Array,
typename T, std::size_t NDims, typename index_t, std::size_t NDimsNew>
multi_array_ref<T, NDimsNew>
reshape(Array<T, NDims>& A, const array<index_t, NDimsNew>& dims) {
multi_array_ref<T, NDimsNew> a(A.origin(), dims);
return a;
}

template <template <typename, std::size_t, typename…> class Array, typename T>
void f(Array<T, 1>& A) {
for (auto it : A) {
// do something with it
std::cout << it <<“”;
}
std::cout << std::endl;
}

template <template <typename, std::size_t, typename…> class Array,
typename T, std::size_t NDims>
void f(Array<T, NDims>& A, long d) {
auto dims = A.shape();
typedef typename std::decay<decltype(*dims)>::type type;

// collapse dimensions [0,d) and (d,Ndims)
array<type, 3> dims3 = {
std::accumulate(dims, dims + d, type(1), std::multiplies<type>()),
dims[d],
std::accumulate(dims + d + 1, dims + NDims, type(1), std::multiplies<type>())
};

// reshape to collapsed dimensions
auto A3 = reshape(A, dims3);

// call f for each slice [i,:,k]
for (auto Ai : A3) {
for (index_t k = 0; k < dims3[2]; ++k) {
auto S = Ai[indices[range()][k]];
f(S);
}
}
}

template <template <typename, std::size_t, typename…> class Array,
typename T, std::size_t NDims>
void f(Array<T, NDims>& A) {
for (long d = NDims; d—; ) {
f(A, d);
}
}

这是测试程序 test.cpp:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25

#include”f.hpp”

int main() {
boost::multi_array<double, 3> A(boost::extents[2][2][3]);
boost::multi_array_ref<double, 1> a(A.data(), boost::extents[A.num_elements()]);
auto Ajk = A[1];
auto Aik = A[boost::indices[range()][1][range()]];

int i = 0;
for (auto& ai : a) ai = i++;

std::cout <<“work on boost::multi_array_ref” << std::endl;
f(a);

std::cout <<“work on boost::multi_array” << std::endl;
f(A);

std::cout <<“work on boost::detail::multi_array:sub_array” << std::endl;
f(Ajk);

std::cout <<“work on boost::detail::multi_array:sub_array” << std::endl;
f(Aik); // wrong result, since reshape() ignores strides!

//f(A[1]); // fails: rvalue A[1] is boost::multi_array_ref<double, 3ul>::reference
}

显然，这种方法存在问题，即当将切片传递给 f() 时，内存不再连续，这会破坏 reshape() 的实现。

似乎更好(更像 C)的方法是从 boost 类型提供的迭代器中构造一个聚合迭代器，因为这会自动处理沿给定维度的非统一步幅。 boost::detail::multi_array::index_gen 看起来很相关，但我不太清楚如何使用它来对维度 d 中的所有切片进行迭代。有什么想法吗？

注意：

SO 上已经有类似的问题，但没有一个让我很满意。我对 N = 3 或 N = 2 的专门解决方案不感兴趣。它必须适用于任何 N。

更新：

这相当于我在 Python 中想要的：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33

def idx_iterator(s, d, idx):
if len(s) == 0:
yield idx
else:
ii = (slice(None),) if d == 0 else xrange(s[0])
for i in ii:
for new_idx in idx_iterator(s[1:], d – 1, idx + [i]):
yield new_idx

def iterator(A, d=0):
for idx in idx_iterator(A.shape, d, []):
yield A[idx]

def f(A):
for d in reversed(xrange(A.ndim)):
for it in iterator(A, d):
print it
print

import numpy as np
A = np.arange(12).reshape((2, 2, 3))

print“Work on flattened array”
f(A.ravel())

print“Work on array”
f(A)

print“Work on contiguous slice”
f(A[1])

print“Work on discontiguous slice”
f(A[:,1,:])

使用 index_gen.hpp 中的功能应该可以以某种方式实现相同的功能，但我仍然无法弄清楚如何。

相关讨论

你看到 MultiArray Concept 下的解释了吗？特别是。 ” 例如，如果 indices 是 index_gen 类型的对象，则以下示例：indices[index_range(0,5)][2][index_range(0,4)]; 具有退化的第二维。从上述规范生成的视图将具有 2 个维度，形状为 5 x 4″
我注意到。请注意我 5 分钟前是如何删除该评论的 :)
@sehe：是的，我读过这个。假设我得到这样一个 5 x 4 视图对象。然后我需要首先迭代 5 维度，沿 4 维度的所有切片应用 f()，反之亦然。在这种情况下，我的方法失败了，因为内存不是连续的。这就是为什么我认为我必须以某种方式利用视图对象提供的迭代器。
您可以使用迭代器迭代视图。如果内存实际上是连续的很重要，那么无论如何您都将不得不复制数据
@sehe：如果工作维度为 d，则将这些迭代器用于 ND 数组会变得复杂。如果内存是连续的，我的方法已经奏效了。如果任何/部分/所有维度不连续，我也希望它能够工作。
你说的很明显。你想要更复杂的东西，所以你不能使用简单的方法。您宁愿编写复杂的迭代，还是使用已实现的库？

我不清楚如何使用”图书馆”。我怀疑 index_gen 是相关的，但我还没有弄清楚如何使用它来实现我想要的。

好的，在花了大量时间研究 boost::multi_array 的实现之后，我现在准备回答我自己的问题：不，在 boost::multi_array 中没有任何规定允许人们沿着任何地方进行迭代，但第一个维度。我能想到的最好方法是构建一个迭代器，手动管理正在迭代的 N-1 索引。这是 slice_iterator.hpp:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71

#include”boost/multi_array.hpp”

template <template <typename, std::size_t, typename…> class Array,
typename T, std::size_t NDims>
struct SliceIterator {
typedef Array<T, NDims> array_type;
typedef typename array_type::size_type size_type;
typedef boost::multi_array_types::index_range range;
typedef boost::detail::multi_array::multi_array_view<T, 1> slice_type;
typedef boost::detail::multi_array::index_gen<NDims, 1> index_gen;

array_type& A;
const size_type* shape;
const long d;
index_gen indices;
bool is_end = false;

SliceIterator(array_type& A, long d) : A(A), shape(A.shape()), d(d) {
int i = 0;
for (; i != d; ++i) indices.ranges_[i] = range(0);
indices.ranges_[i++] = range();
for (; i != NDims; ++i) indices.ranges_[i] = range(0);
}

SliceIterator& operator++() {
// addition with carry, excluding dimension d
int i = NDims – 1;
while (1) {
if (i == d) —i;
if (i < 0) {
is_end = true;
return *this;
}
++indices.ranges_[i].start_;
++indices.ranges_[i].finish_;
if (indices.ranges_[i].start_ < shape[i]) {
break;
} else {
indices.ranges_[i].start_ = 0;
indices.ranges_[i].finish_ = 1;
—i;
}
}
return *this;
}

slice_type operator*() {
return A[indices];
}

// fakes for iterator protocol (actual implementations would be expensive)
bool operator!=(const SliceIterator& r) {
return !is_end;
}

SliceIterator begin() {return *this;}
SliceIterator end() {return *this;}
};

template <template <typename, std::size_t, typename…> class Array,
typename T, std::size_t NDims>
SliceIterator<Array, T, NDims> make_slice_iterator(Array<T, NDims>& A, long d) {
return SliceIterator<Array, T, NDims>(A, d);
}

// overload for rvalue references
template <template <typename, std::size_t, typename…> class Array,
typename T, std::size_t NDims>
SliceIterator<Array, T, NDims> make_slice_iterator(Array<T, NDims>&& A, long d) {
return SliceIterator<Array, T, NDims>(A, d);
}

可以作为

1
2
3

for (auto S : make_slice_iterator(A, d)) {
f(S);
}

适用于我的问题中的所有示例。

我必须说，boost::multi_array 的实现让我非常失望：超过 3700 行的代码应该只是索引管理。特别是只为第一个维度提供的迭代器，与性能实现相差甚远：实际上每一步都进行了多达 3*N + 5 的比较，以确定迭代器是否已经到达末尾(注意我上面的实现通过伪造 operator!=()) 来避免这个问题，这使得这个实现不适合任何东西，除了具有主导最后一维的数组，它可以更有效地处理。此外，该实现没有利用内存中连续的维度。相反，它总是对数组赋值等操作逐维进行，浪费了大量的优化机会。

总之，我发现 numpy 的 N 维数组实现比这个更引人注目。还有 3 个观察结果告诉我 boost::multi_array:

的”放手”

我在网络上的任何地方都找不到 boost::multi_array 的任何严重用例
发展似乎在 2002 年基本停止
StackOverflow 上的这个(和类似的)问题几乎没有引起任何兴趣 ;-)

来源：https://www.codenong.com/23758490/

Iterate over all but d-th dimension of any boost::multi_array

猜你喜欢