Compile Time C++ Tensors via Parameter Packs

The manipulation of tensors (n-dimensional arrays) emerges from several application domains in engineering, computer science and statistics (e.g. discrete CPDs). In C++ this is typically possible with dynamically allocated multidimensional arrays such as the ones from OpenCV (Mat class), but in many occasion there is the interest to rely on the optimizations provided by compile-size arrays.
Two approaches have been proposed: the first is based on the nesting of std::array (here) while the second, used by boost multidimensional library (here),  provides views over single fixed size array. In the latter case the size is fixed at compile-time but the multidim shape is not and there is partial compile-time size verification.

In this post I am addressing full compile-time multidimensional types obtained via C++11 template metaprogramming in particular parameter packs. A class will provide the storage of data as a flat array, and then different views can be generated over this array. In particular the array employs Eigen for efficiently providing mathematical operations.

The core of the library is the representation of multidimensional arrays and their views, based on the element type E and a sequence pairs containing the size of the given dimension, and the step. The pair is represented by a template class called sspair. The offset of a multi-index is obtained by summing the pair-wise products between the step and the index. For row-major ordering each step is the result of the product of higher numbered dimensions.

The resulting type can be used as follows:

MultiDimNRow<double,5,6,7,8> A;
A.limit1<2>(2).setZero(); // fixes dimension 2 to offset 2 and zeros

The code is available on github.

Update 2015-10-15: in the Microsoft GLS (first public on 19th August) the array_view provides multidimensional access to arrays. In particular array_view provides both static and dynamic sizes with support for compact and strided arrays. The multidim discussed here provides only static sizes with the advantage of compile time strided arrays. In addition, due to the mathematical target this solution uses Eigen for the underlying math and some reshape/permute functions.

Update 2015-10-28: look also at similar problem in a talk of CppCon15 "C++ Multi-dimensional Arrays for Computational Physics and Applied Mathematics - Pramod Gupta"

Update 2016-01-29: note that Eigen 3.3 beta contains the unsupported Tensor module developed by Google that provides many interesting features.

Structure

The MultiDimN class is templated over a sequence of sspair types, one per dimension. The number of elements of the array is the product of the steps values of every sspair. Instead of directly expressing the sequence as a parameter pack, the MultiDimN is parametrized against a type_sequence. The motivation is that the type_sequence (abbreviated in TS) can be manipulated in an easier way than the parameter pack, and in particular I am providing the following operations:

• append/prepend a type to a TS
• extract a type from a TS by index
• reverse the order of types
• permute the order of types by indexes
• remove a type by index
• replace a type by index
• remove all types that match a predicate
• perform numeric operations over the TS

Update: The TS can be found also in Boost MPL as TypeList but C++11 allows to make it quite elegant. In comparison to TypeList this implementation provides reverse, permutation by indices and replacement of a type. These operations are indeed specific to the multidimensional array.

Type_sequence interface

template<class...I> struct type_sequence
{ 
  static constexpr int size = sizeof...(I);

  template<class T>
  using append = type_sequence<I..., T>;

  template<class T>
  using prepend = type_sequence<T, I...>;

  template<int looked>
  using pick = typename details::saccessor<looked,0,I...>::type;

  template <class T = void> 
  using inverted = details::invertmake<I...>;

  template <int index> 
  using drop = details::droppermake<index, I...>;

  template <int index, class newtype> 
  using replacetype = details::replacetypemake<index, newtype, I...>;

  template <template <class T> class Pred>
  using removeif = details::removeif_make<Pred,I...>; 

};

Details

The interesting part of the type_sequence is the recursive scheme used for implementing operations such as inverted, drop and removeif. These operations relies on two support types: type_holder concept and std::integral_constant. A type_holder is a type that has to define a type named "type", while the std::integral_constant a value static constexpr of type int named "value".

The following is an example of the drop of type by index, obtained by inheritance recursion pattern and boolean template specialization. The recursion pattern is meant to specify the final output as a type_holder containing the resulting type_sequence, while scanning the indices of the types from 0 to the last one, and, at the same time, bringing down the resulting type_sequence. The template specialization on true/false is used when the current index matches the requested one.

template<bool match, template <class T> class Pred , class OT, class...Is> struct __drophelp;

  /// build a parameter pack with size sizeof(N)-1 removing the j-th
  template <template <class T> class Pred, class OT, class...Is>
  struct dropper;

  /// no more to process
  template <template <class T> class Pred, class...Pre>
  struct dropper<Pred,  type_sequence<Pre...> >
    : type_holder<type_sequence<Pre...> > 
    {
    };

  template <template <class T> class Pred, class...Pre, class x, class...Past>
  struct dropper<Pred,  type_sequence<Pre...>, x, Past...>
    : __drophelp<Pred<x>::value, Pred,type_sequence<Pre...>, x, Past...>
    {
    };

  // positive means remove
  template<template <class T> class Pred,class x, class...Pre, class...Past> 
  struct __drophelp<true,Pred,type_sequence<Pre...>, x, Past...>:
     dropper<Pred,type_sequence<Pre...>,Past... >
  {
  };

  // negative remove
  template<template <class T> class Pred,class x, class...Pre, class...Past> 
  struct __drophelp<false,Pred,type_sequence<Pre...>,x,Past...>:
     dropper<Pred,type_sequence<Pre...,x>, Past... >
  {
  };

template <int j, class...Is>
  using droppermake = typename dropper<j,0,type_sequence<>, Is...>::type;


The final type alias is used to start the recursion scheme with the 0 index and the empty output.

Multidimensional array

The common class is MultiDimNBase and then two classes are exposed, one for the data holder class MultiDimN, based on Eigen, and a view MultiDimNView, that holds simply a pointer to the container. The creation of the MultDimN type is possible in col-major or row-major form, an operation performed by recursion over N... building a type_sequence of sspair. 

template <class T, int...N>
using MultiDimNRow = MultiDimN<T, typename details::rowmajorstepper<N...> >;

The computation of offsets from indexes is provided by nested parameter packs plus an intermediate type that multiplies the step value of sspair by an int:

template<class... TS> struct sspairbyindex {
        template<int ...Args2> struct byindex {
            using type = typename type_sequence<intproductstep<TS,Args2>...>::template sumvalued<int>;
        };
};

template<class... TS,int...I>
    struct offsetcompute<type_sequence<TS...>, I...> {
      using type = typename sspairbyindex<TS...>::template byindex<I...>::type;
    };

A view can be created by fixing the index of a given dimension, or, instead, by reducing the size of a given dimension starting from a given index. The size of the resulting multidim is computed at compile-time while the index can be provided at run-time. The former operation uses drop of the TS to remove the dimension, while the latter uses the replacetype operation updating the given dimension with the new size while preserving the steps.

The squeeze operation, that removes all the singleton dimensions, is provided via the remove_if of the TS. 

Closing

Most relevant operations of multidim has been included, and others can be obtained by direct manipulation of the type_sequence or by similar approaches. Some possible additions:

• Assignment between multidim with same sizes but different steps: boolean predicate on two TS
• Dimensions permutations (aka transpose): simply a reordering of TS, invert provides generalized transposition
• Iteration over all the multdim slices span by the combination of a subset of the indices: recursive enumeration depending on matching indices
• Matrix multiplication extended to multidim: related to the iteration

The dynamically sized version of multidim has been also implemented but not included in the repository: it relies on a shared_ptr of MatrixX of Eigen for storing and sharing data, plus two Eigen integer vectors representing size and steps respectively. The offset computation in this case is the dot product of the indices with the steps vector.

Having both static and dynamic types allows for a single type-definition much like Eigen::Matrix. In addition, due to the increased amount of code generated by the template variants it could be interesting to have a solution in the middle between static and dynamic types: compile-time check of the types and then one single class for the implementation. 

Finally note that the manipulations provided by type_sequence<...> can be applied to other situations and they are not specific of multidim, some of the mentioned operations can be also implemented over integer_sequence. As later discovered check MPL TypeList for a comparison.

Planned Post Updates

Examples of use of the multidim types and their declaration. Example of the accessor via call operator overloading (much like Eigen) taking a std::initializer_list or a parameter pack.

Update on 2015-10-19

• added permutation type-sequence for supporting dimension permutation
• added reshape (col-wise or row-wise) returned as view