## 13 January 2010

### Providing fill and for_each algorithms for Boost.MultiArray

I sometimes want to fill a Boost.MultiArray with a particular value. For concrete implementations like boost::multi_array or boost::multi_array_ref, the task isn't too hard and std::fill works just fine. But getting an arbitrary-dimensional version working that only uses the MultiArray concept takes a little thought. This is important when working with views. I additionally want the solution to take advantage of the contiguous storage guarantees made by multi_array and multi_array_ref.

Because fill looks like for_each with an assignment-like UnaryFunction, I decided to write a MultiArray-specific for_each implementation. First comes a helper functor recursively templated on the dimensionality:

namespace { // anonymous

template<std::size_t D>
struct for_each_functor {

BOOST_STATIC_ASSERT(D != 0); // Nonsensical behavior for zero dimensions
BOOST_STATIC_ASSERT(D > 1);  // Ensure not instantiated for specialized values

// for details about why MultiArray::iterator::reference is used below.
template<class MultiArray,
class UnaryFunction >
void operator()(MultiArray &x, UnaryFunction &f) const {
for_each_functor<D-1> functor;
for (typename MultiArray::iterator i = x.begin(); i != x.end(); ++i) {
typename MultiArray::iterator::reference ri = *i;
functor(ri,f);
}
}
};

template<> struct for_each_functor<1> {
template<class MultiArray,
class UnaryFunction >
void operator()(MultiArray &x, UnaryFunction &f) const {
std::for_each(x.begin(), x.end(), f);
}
};

} // namespace anonymous

The functor uses std::for_each for the one dimensional base case. The BOOST_STATIC_ASSERTs are there just to keep me from shooting myself in the foot. A smarter version might specialize for the two- and three-dimensional cases, check the underlying strides, and then be certain to iterate in the faster directions first.

Now that we've got for_each_functor that works on the general MultiArray concept, it's straightforward to implement for_each:

/**
* Invoke \c f for each element in MultiArray \c x.  The
* order in which the invocations take place is undefined.
*
* @param x MultiArray over which to iterate.
* @param f UnaryFunction to invoke on each element of \c x.
*
* @return \c f.
*
* @see SGI's <a href="http://www.sgi.com/tech/stl/UnaryFunction.html">
*/
template<class MultiArray,
class UnaryFunction>
inline
UnaryFunction for_each(MultiArray &x,
UnaryFunction f) {
for_each_functor<MultiArray::dimensionality>()(x,f);
return f;
}

Please note that my for_each type signature is not identical with the STL algorithm's signature.

Since we know that multi_array and multi_array_ref have contiguous storage, we can specialize our version of for_each to linearly walk their underlying memory using the usual algorithms.

/**
* Invoke \c f for each element in <tt>boost::multi_array</tt> \c x.  The order
* in which the invocations take place is undefined.  This specialization
* takes advantage of <tt>boost::multi_array</tt>'s contiguous storage.
*
* @param x <tt>boost::multi_array</tt> over which to iterate.
* @param f UnaryFunction to invoke on each element of \c x.
*
* @return \c f.
*
* @see SGI's <a href="http://www.sgi.com/tech/stl/UnaryFunction.html">
*/
template<class ValueType, std::size_t NumDims, class Allocator,
class UnaryFunction>
inline
UnaryFunction for_each(boost::multi_array<ValueType,NumDims,Allocator> &x,
UnaryFunction f) {
return std::for_each(x.data(), x.data() + x.num_elements(), f);
}

/**
* Invoke \c f for each element in <tt>boost::multi_array_ref</tt> \c x.  The
* order in which the invocations take place is undefined.  This specialization
* takes advantage of <tt>boost::multi_array_ref</tt>'s contiguous storage.
*
* @param x <tt>boost::multi_array_ref</tt> over which to iterate.
* @param f UnaryFunction to invoke on each element of \c x.
*
* @return \c f.
*
* @see SGI's <a href="http://www.sgi.com/tech/stl/UnaryFunction.html">
*/
template<class ValueType, std::size_t NumDims,
class UnaryFunction>
inline
UnaryFunction for_each(boost::multi_array_ref<ValueType,NumDims> &x,
UnaryFunction f)
{
return std::for_each(x.data(), x.data() + x.num_elements(), f);
}

With for_each complete, we can move on to implementing fill by mapping an assignment-like functor over over a MultiArray.

In my particular application I have MultiArrays with complex-valued elements from a couple of different numeric libraries (e.g. std::complex<T>, T[2]) and I'd like to be able to fill a MultiArray with a value for which a member operator= may not be available. To accommodate this need, I created a templated assignment functor which can be specialized using schmancy boost::enable_if techniques based on both the target and source of the assignment:

/**
* A functor that performs assignment to type \c Target from type \c Source.
* \c Target must be assignable from \c Source.
*
* The \c Enable template parameter allows using <tt>boost::enable_if</tt> to
* specialize or partially specialize the functor per <tt>enable_if</tt>'s
* documentation section <a
* href="http://www.boost.org/doc/libs/1_41_0/libs/utility/enable_if.html">
* 3.1: Enabling template class specializations</a>.
*/
template<class Target, class Source, class Enable = void>
struct assign {

/**
* Create an instance which assigns \c s when applied.
*
* @param s source of assignment operation occurring via
*          <tt>operator()</tt>.
*/
assign(const Source &s) : s_(s) {}

/**
* Assign the value provided at construction to \c t.
*
* @param t to be assigned.
*/
void operator()(Target& t) const { t = s_; }

private:
const Source &s_; /**< Source for assignment operations */
};

The default implementation above makes the extra parameters seem overkill. I'll give an example where it is useful later.

Using the assignment functor and for_each, here's our one-line fill:

/**
* Fill MultiArray \c x with the value \c v.  MultiArray <tt>x</tt>'s elements
* must be assignable from \c v.  The underlying assignment uses
* ::assign to allow specializations of that functor to be found.
*
* @param x MultiArray to fill.
* @param v Value with which to fill \c x.
*/
template<class MultiArray, class V>
void fill(MultiArray &x, const V &v) {
for_each(x, assign<typename MultiArray::element,V>(v));
}

Finally, to show the spurious-looking assign functor template parameters are useful, here's an assign functor for the complex-valued use case I mentioned above:

/**
* A specialization of the assign functor to properly handle the case where \c
* Target is a recognized complex type according to
* ::suzerain::complex::traits::is_complex.  It uses
* ::suzerain::complex::assign_complex to perform the assignment, and therefore
* supports all types that \c assign_complex does.
*/
template<class Target, class Source>
struct assign<
Target,
Source,
typename boost::enable_if<
::suzerain::complex::traits::is_complex<Target>
>::type >
{
/**
* Create an instance which assigns \c s when applied.
*
* @param s source of assignment operation occurring via
*          <tt>operator()</tt>.
*/
assign(const Source &s) : s_(s) {};

/**
* Assign the value provided at construction to \c t.
*
* @param t to be assigned.
*/
void operator()(Target& t) const {
::suzerain::complex::assign_complex(t, s_);
}

private:
const Source &s_; /**< Source for assignment operations */
};

That very last code snippet is missing details and won't compile. Note that the suzerain::complex::assign_complex invocation may itself be a template and it receives appropriate type information. This turned out to be important when I wanted to fill a complex-valued MultiArray with a real-valued scalar. I needed two versions of assign_complex: one that took as a Source a recognized complex type and one that took a scalar type. boost::enable_if and boost::disable_if made it easy to provide an assign_complex template that did the right thing.

Given the right #includes in the right places, the other snippets should all compile. Please let me know if they do not.

This post topic came from a question I asked on boost-users a couple of days back. Hopefully someone else finds this post and can use this fill routine. Thank you to Ronald Garcia for his response and for maintaining Boost.MultiArray.