How to Traverse a Boost::Multi_Array

Iterate through a boost::multi_array view

Here is one way to do this:

#include <iostream>
#include <boost/multi_array.hpp>

// Functor to iterate over a Boost MultiArray concept instance.
template<typename T, typename F, size_t Dimensions = T::dimensionality>
struct IterateHelper {
   void operator()(T& array, const F& f) const {
      for (auto element : array)
         IterateHelper<decltype(element), F>()(element, f);
   }
};

// Functor specialization for the final dimension.
template<typename T, typename F>
struct IterateHelper<T, F, 1> {
   void operator()(T& array, const F& f) const {
      for (auto& element : array)
         f(element);
   }
};

// Utility function to apply a function to each element of a Boost
// MultiArray concept instance (which includes views).
template<typename T, typename F>
static void iterate(T& array, const F& f) {
   IterateHelper<T, F>()(array, f);
}

int main() {
   boost::multi_array<char, 2> a{boost::extents[2][6]};

   a[0][0] = 'B';
   a[0][1] = 'o';
   a[0][2] = 'o';
   a[0][3] = 's';
   a[0][4] = 't';
   a[0][5] = '\0';

   a[1][0] = 'L';
   a[1][1] = 'i';
   a[1][2] = 'o';
   a[1][3] = 'n';
   a[1][4] = 's';
   a[1][5] = '\0';

   typedef boost::multi_array<char, 2>::array_view<2>::type array_view;
   typedef boost::multi_array_types::index_range range;
   array_view b = a[boost::indices[range{0,2}][range{0,4}] ];

   // Use the utility to apply a function to each element.
   iterate(b, [](char& c) {
      std::cout << c << std::endl;
   });

   return 0;
};

The code above defines a utility function iterate(), to which you pass an object satisfying the Boost MultiArray concept (which includes views) and a function to apply to each element. The utility function works by using a Functor that iterates over each dimension recursively.

CoLiRu

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

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

Dimension-independent loop over boost::multi_array?

Ok, based on the Google groups discussion already mentioned in one of the comments and on one of the examples from the library itself, here is a possible solution that lets you iterate over all values in the multi-array in a single loop and offers a way to retrieve the index for each of these elements (in case this is needed for some other stuff, as in my scenario).

#include <iostream>
#include <boost/multi_array.hpp>
#include <boost/array.hpp>

const unsigned short int DIM = 3;
typedef double tValue;
typedef boost::multi_array<tValue,DIM> tArray;
typedef tArray::index tIndex;
typedef boost::array<tIndex, DIM> tIndexArray;

tIndex getIndex(const tArray& m, const tValue* requestedElement, const unsigned short int direction)
{
  int offset = requestedElement - m.origin();
  return(offset / m.strides()[direction] % m.shape()[direction] +  m.index_bases()[direction]); 
}

tIndexArray getIndexArray( const tArray& m, const tValue* requestedElement )
{
  tIndexArray _index;
  for ( unsigned int dir = 0; dir < DIM; dir++ )
  {
    _index[dir] = getIndex( m, requestedElement, dir );
  }

  return _index;
}

int main()
{ 
  double* exampleData = new double[24];
  for ( int i = 0; i < 24; i++ ) { exampleData[i] = i; }

  tArray A( boost::extents[2][3][4] );
  A.assign(exampleData,exampleData+24);

  tValue* p = A.data();
  tIndexArray index;
  for ( int i = 0; i < A.num_elements(); i++ )
  {
    index = getIndexArray( A, p );
    std::cout << index[0] << " " << index[1] << " " << index[2] << " value = " << A(index) << "  check = " << *p << std::endl;
    ++p;
  }

  return 0;
}

The output should be

0 0 0 value = 0 check = 0
0 0 1 value = 1 check = 1
0 0 2 value = 2 check = 2
0 0 3 value = 3 check = 3
0 1 0 value = 4 check = 4
0 1 1 value = 5 check = 5
0 1 2 value = 6 check = 6
0 1 3 value = 7 check = 7
0 2 0 value = 8 check = 8
0 2 1 value = 9 check = 9
0 2 2 value = 10 check = 10
0 2 3 value = 11 check = 11
1 0 0 value = 12 check = 12
1 0 1 value = 13 check = 13
1 0 2 value = 14 check = 14
1 0 3 value = 15 check = 15
1 1 0 value = 16 check = 16
1 1 1 value = 17 check = 17
1 1 2 value = 18 check = 18
1 1 3 value = 19 check = 19
1 2 0 value = 20 check = 20
1 2 1 value = 21 check = 21
1 2 2 value = 22 check = 22
1 2 3 value = 23 check = 23

so the memory layout goes from the outer to the inner indices. Note that the getIndex function relies on the default memory layout provided by boost::multi_array. In case the array base or the storage ordering are changed, this would have to be adjusted.

Iterating over the dimensions of a boost::multi_array

Normally you can do it with boost::multi_array, heres some sample code to create a 2D view of a 3D multi_array:

#include "boost/multi_array.hpp"
#include <cassert>
#include <iostream>

int main()
{
    typedef boost::multi_array<double, 3> array_type;
    typedef array_type::index index;

    array_type myArray3D(boost::extents[3][4][2]);

    // Assign values to the elements
    int values = 0;
    for (index i = 0; i != 3; ++i)
    {
            for (index j = 0; j != 4; ++j)
            {
                    for (index k = 0; k != 2; ++k)
                    {
                            myArray3D[i][j][k] = values++;
                    }
            }
    }
    // Verify values
    int verify = 0;
    for (index i = 0; i != 3; ++i)
    {
            for (index j = 0; j != 4; ++j)
            {
                    for (index k = 0; k != 2; ++k)
                    {
                            std::cout << "[" << i << "]";
                            std::cout << "[" << j << "]";
                            std::cout << "[" << k << "] = ";
                            std::cout << myArray3D[i][j][k] << std::endl;
                            assert(myArray3D[i][j][k] == verify++);
                    }
            }
    }

    typedef boost::multi_array_types::index_range range;
    array_type::index_gen indices;

    // Create a new view with 2 dimentions fixing the 2nd dimention to 1
    array_type::array_view<2>::type myView =
    myArray3D[indices[range()][1][range()]];
    std::size_t numDims = myView.size();
    std::cout << "numDims = " << numDims << std::endl;

    for (index i = 0; i != 3; ++i)
    {
            for (index j = 0; j != 2; ++j)
            {
                    std::cout << "[" << i << "]";
                    std::cout << "[" << j << "] = ";
                    std::cout << myView[i][j] << std::endl;
            }
    }

    return 0;
}

Boost::multi_array looping

The index_bases member function returns a container with each dimension's index base. The shape member function returns a container with each dimension's extent (size). You can use both of these to determine the range of indices for each dimension:

typedef boost::multi_array<int, 3> array_type;

void printArray(const array_type& a)
{
    // Query extents of each array dimension
    index iMin = a.index_bases()[0];
    index iMax = iMin + a.shape()[0] - 1;
    index jMin = a.index_bases()[1];
    index jMax = jMin + a.shape()[1] - 1;
    index kMin = a.index_bases()[2];
    index kMax = kMin + a.shape()[2] - 1;

    for (index i=iMin; i<=iMax; ++i)
    {
        for (index j=jMin; j<=jMax; ++j)
        {
            for (index k=kMin; k<=kMax; ++k)
            {
                std::cout << a[i][j][k] << " ";
            }
        }
    }
}

int main()
{
    typedef array_type::extent_range range;
    typedef array_type::index index;
    array_type::extent_gen extents;

    // Extents are hard-coded here, but can come from user/disk.
    array_type a(extents[2][range(1,4)][range(-1,3)]);

    // Populate array with values...

    // Pass the array to a function. The function will query
    // the extents of the given array.
    print(a);

    return 0;
}

Boost C++ - dynamically iterating over multi-array

I looked into the boost code for boost::multi_array_types::index_range, and found that this is not possible. The class has only three members to store index values - start, finish, and stride. It can't store a more complex set of values.

Since the number of columns that I need is dynamic, I used a vector of sub-arrays (array_view),

vector<boost::multi_array_ref<double, 2>::array_view<2>::type::const_reference

and just added to the vector as needed.

---

Related Topics