Implement STL functions in variadic template

Question

I have been working on a small project to get up to speed with variadic templates. I implemented a small multidimensional array. I would now like to define a function that operates on the nearest neighbours of a given position -- is there an elegant way to retrieve the values of the neighbours of a given position in my array?

template<class T, size_t size, size_t... sizes>
struct MArr {
    typedef std::array<typename MArr<T, sizes...>::type, size> type;

    std::array<MArr<T, sizes...>,size> data;

    MArr<T, sizes...>& operator[](int i) {
        return data[i];
    }

};

template<class T, size_t size>
struct MArr<T, size> {
    typedef std::array<T, size> type;
    type data;

    T& operator[](int i) {
       return data[i];
    }

};

Addendum: It is somewhat clear to me how to loop all elements by employing recursion, e.g. applying an arbitrary function to a hyperdim. array of ints:

template <typename T, size_t size>
void func(MArr<T, size>& arr, std::function<void(int &)> f) {
    for (int i = 0; i < size; i++) {
       f(arr[i]);
    }
}


template <typename T, size_t size0, size_t size1, size_t ...sizes>
void func(MArr<T, size0, size1, sizes...>& arr, std::function<void(int &)> f) {
    for (int i =0; i < size0; i++) {
        func<T, size1, sizes...>(arr[i], f);
    }
}

I am curious how I would generate something like func(arr[i+1,j,k,…],arr[i-1,j,k,…],arr[i,j+1,k,…],arr[i,j-1,k,…],…) and a given func (which, say adds the respective elements). As I said, quite new to variadic templates and I have the feeling I don't have the right mindset yet …

Answer 1

You may do something like (code use folding expression from C++17, but can be written in C++11):

template <std::size_t I, typename F, std::size_t ... Is, typename Tuple>
void helper(F&& f, std::index_sequence<Is...>, const Tuple& t)
{
    f((std::get<Is>(t) - (Is == I))...);
    f((std::get<Is>(t) + (Is == I))...);
}

template <typename F, std::size_t ... Is, typename Tuple>
void helper(F&& f, std::index_sequence<Is...> Seq, const Tuple& t)
{
    (helper<Is>(std::forward<F>(f), Seq, t), ...);
}

template <typename F, typename ... Ts>
void apply_to_neighboor(F&& f, Ts... indexes)
{
    helper(std::forward<F>(f), std::index_sequence_for<Ts...>(), std::tie(indexes...));
}

Demo

If you want to retrieve all neighbours, you might change above code to:

template <std::size_t I, std::size_t ... Is, typename Tuple>
auto helper1(std::index_sequence<Is...>, const Tuple& t)
{
    return std::make_pair(std::make_tuple((std::get<Is>(t) - (Is == I))...),
                          std::make_tuple((std::get<Is>(t) + (Is == I))...));
}

template <std::size_t ... Is, typename Tuple>
auto helper(std::index_sequence<Is...> Seq, const Tuple& t)
{
    return std::tuple_cat(helper1<Is>(Seq, t)...);
}

template <typename F, typename ... Ts>
void apply_to_neighboor(F&& f, Ts... indexes)
{
    std::apply(std::forward<F>(f),
               helper(std::index_sequence_for<Ts...>(), std::tie(indexes...)));
}

Demo

Answer 2

Here is a possible implementation that rely on your MArr providing a custom subscript operator:

// Subscript using an array of size_t
T operator[](const size_t (&idx)[sizeof... (sizes) + 1]) const;

The idea is as follows:

1. You generate a sequence from 0 to 2 * D * D - 1 (where D is the number of dimensions).
2. Using this sequence, you generate all the indices in a single 1D array x + 1, y, z, x - 1, y, z, ..., x, y, z + 1, x, y, z - 1.
3. You treat this huge array as an array of 2 * D arrays of D values (size_t [2 * D][D]).
4. You use this array to index your multi-dimensional array using the newly defined subscript operator.

And here would be the implementation:

template <class T, std::size_t... Sizes, class F, class... Is,
          std::size_t... IsP2N, std::size_t... IsP2NN>
auto around_impl(MArr<T, Sizes...> const& arr,
                 std::index_sequence<IsP2N...>,
                 std::index_sequence<IsP2NN...>,
                 F &&f, Is&&... pt) {
    const std::size_t pts[] = {(std::size_t)pt... };
    const std::size_t pts2[2 * sizeof...(Sizes)][sizeof...(Sizes)] = {
        (pts[IsP2NN % sizeof...(Sizes)] 
         + (1 - 2 * ((IsP2NN / sizeof...(Sizes)) % 2)) 
            * (IsP2NN % sizeof...(Sizes) == IsP2NN / (2 * sizeof...(Sizes))))...

    };
    return f(arr[pts2[IsP2N]]... );
}

template <class T, std::size_t... Sizes, class F, class... Is>
auto around(MArr<T, Sizes...> const& arr, F &&f, Is&&... pt) {
    return around_impl(arr,
                       std::make_index_sequence<2 * sizeof...(Sizes)>{},
                       std::make_index_sequence<2 * sizeof...(Sizes) * sizeof...(Sizes)>{},
                       std::forward<F>(f),
                       std::forward<Is>(pt)... );
}

Small notes on the implementation:

• The computation that takes place inside the brackets is explained at the end of this answer.
• The initialization of a 2D array from a single brace-enclosed list is valid C++ (as far as I know).
• Most of this is optimized out by clang with -O3.

You can then call:

// 3D array
MArr<int, 4, 4, 4> arr;

// Function to call:
auto myFun = [](int a, int b, int c, int d, int e, int f) { 
    return a + b + c + d + e + f;
};

// Call around:
around(arr, myFun, x, y, z);

To call:

myFun(arr[x + 1, y, z], arr[x - 1, y, z], 
      arr[x, y + 1, z], arr[x, y - 1, z],
      arr[x, y, z + 1], arr[x, y, z - 1]);

---

Explanation of:

(pts[IsP2NN % sizeof...(Sizes)] 
 + (1 - 2 * ((IsP2NN / sizeof...(Sizes)) % 2)) 
    * (IsP2NN % sizeof...(Sizes) == IsP2NN / (2 * sizeof...(Sizes))))...

• IsP2NN goes from 0 to 2 * D * D - 1.
• pts[...] corresponds to x, y, z, x, y, z, ..., and so on.
• 1 - ... parts generates 1, 1, 1, -1, -1, -1, 1, 1, 1, -1, -1, -1, and so on... It alternates between sequences of 1 and -1.
• * (...) will nullify the previous sequence when the dimension does not match, it generates 1, 0, 0, 1, 0, 0, 0, 1, 0, 0, 1, 0, and so on...

Here is a "matrix" view of what happens:

Is          0   1   2   3   4   5   6   7   8   9  10  11  12  13  14  15  16  17

Is % 3      0   1   2   0   1   2   0   1   2   0   1   2   0   1   2   0   1   2
Is / 3 % 2  0   0   0   1   1   1   0   0   0   1   1   1   0   0   0   1   1   1
Is / 6      0   0   0   0   0   0   1   1   1   1   1   1   2   2   2   2   2   2

pts  --     x   y   z   x   y   z   x   y   z   x   y   z   x   y   z   x   y   z
     --     1   0   0  -1   0   0   0   1   0   0  -1   0   0   0   1   0   0  -1
     --     1   0   0   1   0   0   0   1   0   0   1   0   0   0   1   0   0   1
---------------------------------------------------------------------------------
          x+1   y   z x-1   y   z   x y+1   z   x y-1   z   x   y z+1   x   y z-1

Answer 3

First I'd create an index type.

template<std::size_t N>
using index=std::array<std::size_t, N>;

template<class T, size_t size, size_t... sizes>
struct MArr {
  using my_index=index<sizeof...(sizes)+1>;

and an index ref. Either use gsl::span or write your own:

namespace utility {
  template<class It>
  struct range {
    It b,e;
    It begin()const{return b;}
    It end()const{return e;}
    std::size_t size()const{return std::distance(begin(),end());}
    bool empty()const{return begin()==end();}

    using reference=typename std::iterator_traits<It>::reference;

    reference front()const{return *begin();}
    reference back()const{return *std::prev(end());}
  };
  template<class T>
  struct span:range<T*> {
    span(T* s, T*f):range<T*>{s,f}{}
    span():span(nullptr,nullptr){}
    span(T*s,std::size_t len):span(s,s+len){}
    T&operator[](std::size_t i)const{return this->begin()[i];}
    T* data()const{return this->begin();}
    span without_front(std::size_t n=1)const{ return {this->begin()+(std::min)(n,this->size()), end()}; }
    span without_back(std::size_t n=1)const{ return {this->begin(), end()-(std::min)(n,this->size())}; }
    span only_front(std::size_t n=1)const{ return {this->begin(),this->begin()+(std::min)(n,this->size())}; }
    span only_back(std::size_t n=1)const{ return {end()-(std::min)(n,this->size()),end()}; }
    span mid(std::size_t start, std::size_t len)const{ return without_front(start).only_front(len); }

    template< class U >
    using compatible=std::integral_constant<bool, std::is_convertible<U*,T*>{}&&(sizeof(U)==sizeof(T))>;
    template<class R>
    using compatible_range=compatible< std::decay_t<decltype( *std::declval<R>().data() )> >;
    template<class C,
      std::enable_if_t< compatible_range< C& >, bool> =true, // has .data() that returns good type
      std::enable_if_t< !std::is_same<span, std::decay_t<C>>{}, bool> =true // not own type
    >
    span(C&& c): span( c.data(), c.size() ){}
    template<std::size_t N>
    span( T(&arr)[N] ):span(arr, N){}
    // etc
  };
}

once you have span (either like above, or from gsl), the code becomes much cleaner:

template<std::size_t N>
using index=std::array<std::size_t, N>;

using index_cref=utility::span<std::size_t const>;
using index_ref=utility::span<std::size_t>;

template<class T, size_t size, size_t... sizes>
struct MArr {
  using my_index=index<sizeof...(sizes)+1>;
  T& operator[]( index_cref r ){ return data[r.front()][ r.without_front() ]; }

and in the one-dim case

T& operator[](index_cref r) {
   return data[r.front()];
}

once we have done this your problem becomes easier.

First we rewrite your func to iterate over values of an index. This can be done the way you did it, where you pass in a callback, or you could implement next which advances an index.

bool next_index( index_ref r, index_cref bounds ){
  if (r.empty()||bounds.empty()) return false;
  ++r.back();
  if (r.back()!=bounds.back()) return true;
  r.back()=0;
  return next_index( r.without_back(), bounds.without_back() );
}

now iterating can look like this:

template<class MArr, class F>
void foreach_index( MArr const&, F&& f ){
  using index=typename MArr::index;
  index const bounds = MArr::bounds(); // todo: write
  index cur = {{0}};
  do {
    f(cur);
  } while( next_index(cur, bounds) );
}

which is cleaner, simpler and more efficient than your version (no type erasure).

Nearest neighbour can be written in terms of index_ref easily.

template<class F>
void foreach_neighbour( index_ref where, index_cref bounds, F&& f ){
  for(std::size_t i=0; i<(std::min)(where.size(),bounds.size());++i){
    if (where[i]){ where[i]--; f(where); where[i]++; }
    if (where[i]+1<bounds[i]) { where[i]++; f(where); where[i]--; }
  }
}