很多时候,人们想要f()
沿d
一N
维数组的维度应用操作A
。这意味着遍历所有剩余的维度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
,理想情况下也适用于 rvalue 类型,例如boost::multi_array_ref<T, NDims>::reference
.
我能想到的最好的方法是实现一个reshape()
函数,该函数可用于将 ND 数组重塑为 3D 数组,这样工作维度始终是中间维度。这里是f.hpp
:
#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
:
#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
or的专门解决方案不感兴趣N = 2
。它必须适用于任何N
.
更新:
这是我在 Python 中想要的等价物:
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
,但我仍然无法弄清楚如何。