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

    // See http://groups.google.com/group/boost-list/browse_thread/thread/e16f32c4411dea08
    // 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">
 *      UnaryFunction</a> concept for more information.
 */
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">
 *      UnaryFunction</a> concept for more information.
 */
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">
 *      UnaryFunction</a> concept for more information.
 */
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.