Итерировать по всем измерениям, кроме d-го, любого увеличения :: multi_array

Довольно часто хочется применить операцию 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в идеале также для таких типов значений, как boost::multi_array_ref<T, NDims>::reference,

Лучшее, что я мог придумать, — это реализация reshape() функция, которая может использоваться для преобразования массива ND в трехмерный массив, так что рабочий размер всегда является средним. Вот 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::detail::multi_array::index_gen выглядит актуально, но мне не совсем понятно, как это можно использовать для создания итератора по всем слайсам в измерении d, Есть идеи?

Замечания:

Уже есть похожие вопросы по SO, но ни один из них не был достаточно удовлетворительным для меня. Я не заинтересован в специализированных решениях для N = 3 или же 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, но я до сих пор не смог понять, как.