Iterate over all but d-th dimension of any boost::multi_array
Question: Quite often one wants to apply operation `f()` along dimension `d` of an
`N`-dimensional array `A`. This implies looping over all remaining dimensions
of `A`. I tried to figure out if `boost::multi_array` was capable of this.
Function `f(A)` should work on all varieties of `boost::multi_array`,
including `boost:multi_array_ref`, `boost::detail::multi_array::sub_array`,
and `boost::detail::multi_array::array_view`, ideally also for the rvalue
types such as `boost::multi_array_ref<T, NDims>::reference`.

The best I could come up with is an implementation of a `reshape()` function
that can be used to reshape the ND array into a 3D array, such that the
working dimension is always the middle one. Here is `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);
        }
    }
    

This is the test program `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
    }
    

Clearly, there are problems with this approach, namely when a slice is passed
to `f()`, such that the memory is no longer contiguous, which defeats the
implementation of `reshape()`.

It appears a better (more C++-like) way would be to construct an aggregate
iterator out of the iterators that the boost types provide, since this would
automatically take care of non-unity strides along a given dimension.
`boost::detail::multi_array::index_gen` looks relevant, but it is not quite
clear to me how this can be used to make an iterator over all slices in
dimension `d`. Any ideas?

_Note:_

There are similar questions already on SO, but none was quite satisfactory to
me. I am not interested in specialized solutions for `N = 3` or `N = 2`. It's
got to work for any `N`.

_Update:_

Here is the equivalent of what I want in 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,:])
    

The same should somehow be possible using the functionality in
`index_gen.hpp`, but I have still not been able to figure out how.


Answer: Ok, after spending a significant amount of time studying the implementation of
`boost::multi_array`, I am now ready to answer my own question: No, there are
no provisions anywhere in `boost::multi_array` that would allow one to iterate
along any but the first dimension. The best I could come up with is to
construct an iterator that manually manages the `N-1` indices that are being
iterated over. Here is `slice_iterator.hpp`:

    
    
    #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);
    }
    

It can be used as

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

and works for all examples in my question.

I must say that `boost::multi_array`'s implementation was quite disappointing
to me: Over 3700 lines of code for what should be little more than a bit of
index housekeeping. In particular the iterators, which are only provided for
the first dimension, aren't anywhere near a performance implementation: There
are actually up to `3*N + 5` comparisons carried out at each step to decide
whether the iterator has arrived at the end yet (note that my implementation
above avoids this problem by faking `operator!=()`), which makes this
implementation unsuitable for anything but arrays with a dominant last
dimension, which is handled more efficiently. Moreover, the implementation
doesn't take advantage of dimensions that are contiguous in memory. Instead,
it always proceeds dimension-by-dimension for operations such as array
assignment, wasting significant optimization opportunities.

In summary, I find `numpy`'s implementation of an N-dimensional array much
more compelling than this one. There are 3 more observations that tell me
"Hands Off" of `boost::multi_array`:

  * I couldn't find any serious use cases for `boost::multi_array` anywhere on the web
  * Development appears to have essentially stopped in 2002
  * This (and similar) questions on StackOverflow have hardly generated any interest ;-)

