TRIOT: iterating over (variadic !) tensors without branching
TRIOT stands for Template Recursive Iteration Over Tensors and promises to iterate over tensors of arbitrary dimensions without branching.

Note: TRIOT as-is does not support strided iterations and must adapted for that (can it?)

# Sources
[Web page](https://alg.cs.umt.edu/triot.html)
[Stack Overflow](https://stackoverflow.com/questions/43356889/a-fast-way-to-iterate-through-multidimensional-arrays-with-different-shapes)
Paper: TRIOT: Faster tensor manipulation in C++11 [1608.00099](https://arxiv.org/pdf/1608.00099)
Git [repo](https://bitbucket.org/orserang/triot/)

# Benchmarks (Arraymancer is similar to tuple iteration)

Note: Numpy benchmark was using fixed C for loops to not penalize Python interpreter (hence should be ignored).

![image](https://user-images.githubusercontent.com/22738317/33306417-cd047e26-d412-11e7-9d87-65feb9c51afb.png)
![image](https://user-images.githubusercontent.com/22738317/33306422-d458f0da-d412-11e7-860f-ddad445f97f8.png)
![image](https://user-images.githubusercontent.com/22738317/33306455-f572a43c-d412-11e7-9bf3-8bb0271205c4.png)

# Implementation
Note, from the StackOverflow post, it seems like it can also support a variadic number of tensors


Unroll shape at compile-time
```C++
template <unsigned int DIMENSION>
inline unsigned long tuple_to_index_fixed_dimension(const_tup_t tup, const_tup_t shape) {
  unsigned long res = 0; unsigned int k;
  for (k=0; k<DIMENSION−1; ++k) {
    res += tup[k];
    res *= shape[k+1]; 
  }
  res += tup[k];
  return res; 
}

template <unsigned int DIMENSION, unsigned int CURRENT> 
class ForEachFixedDimensionHelper {
public:
  template <typename FUNCTION, typename ...TENSORS>
  inline static void apply(tup_t counter, const_tup_t shape, FUNCTION function, TENSORS & ...args) {
    for (counter[CURRENT]=0; counter[CURRENT]<shape[CURRENT]; ++counter[CURRENT]) 
      ForEachFixedDimensionHelper<DIMENSION−1, CURRENT+1>::template apply<FUNCTION, TENSORS...>(counter, shape, function, args...);
  } 
};

template <unsigned int CURRENT>
class ForEachFixedDimensionHelper<1u, CURRENT> { 
public:
  template <typename FUNCTION, typename ...TENSORS>
  inline static void apply(tup_t counter, const_tup_t shape, FUNCTION function, TENSORS & ...args) {
    for (counter[CURRENT]=0; counter[CURRENT]<shape[CURRENT]; ++counter[CURRENT]) 
      function(args[tuple_to_index_fixed_dimension<CURRENT+1>(counter, args.data_shape())]...); /* tensor.data_shape() is an accessor for returning the shape member. */
  } 
};

template <unsigned char DIMENSION> 
class ForEachFixedDimension { 
public:
  template <typename FUNCTION, typename ...TENSORS>
  inline static void apply(const_tup_t shape, FUNCTION function, TENSORS & ...args) {
    unsigned long counter[DIMENSION];
    memset(counter, 0, DIMENSION*sizeof(unsigned long)); 
    ForEachFixedDimensionHelper<DIMENSION,0>::template apply<FUNCTION, TENSORS...>(counter, shape, function, args...);
  } 
};
```

Look-up at runtime
```C++
typedef unsigned int TEMPLATE_SEARCH_INT_TYPE;
template <TEMPLATE_SEARCH_INT_TYPE MINIMUM, TEMPLATE_SEARCH_INT_TYPE MAXIMUM, template <TEMPLATE_SEARCH_INT_TYPE> class WORKER> 
class LinearTemplateSearch {
public:
  template <typename...ARG_TYPES>
  inline static void apply(TEMPLATE_SEARCH_INT_TYPE v, ARG_TYPES && ... args) {
    if (v == MINIMUM)
      WORKER<MINIMUM>::apply(std::forward<ARG_TYPES>(args)...);
    else
      LinearTemplateSearch<MINIMUM+1, MAXIMUM, WORKER>::apply(v, std::forward<ARG_TYPES>(args)...);
  } 
};

template <TEMPLATE_SEARCH_INT_TYPE MAXIMUM, template <TEMPLATE_SEARCH_INT_TYPE> class WORKER >
class LinearTemplateSearch<MAXIMUM, MAXIMUM, WORKER> { 
public:
  template <typename...ARG_TYPES>
  inline static void apply(TEMPLATE_SEARCH_INT_TYPE v, ARG_TYPES && ... args) {
    assert(v == MAXIMUM);
    WORKER<MAXIMUM>::apply(std::forward<ARG_TYPES>(args)...); 
  }
};
```

# Questions
There seems to be 2 constraints:
- DIMENSIONS, which seems similar to Arraymancer MAXRANK
- MINIMUM / MAXIMUM which would be 0 .. 2^31 - 1. The space consumed for precomputation might produce a huge binary. The linear search at compile-time might be prohibitive for fast development.

# Arraymancer implementation:

## Reference implementation
```Nim
template strided_iteration[T](t: Tensor[T], strider: IterKind): untyped =
  ## Iterate over a Tensor, displaying data as in C order, whatever the strides.

  ## Iterator init
  var coord = newSeq[int](t.rank) # Coordinates in the n-dimentional space
  var backstrides: seq[int] = @[] # Offset between end of dimension and beginning
  for i,j in zip(t.strides,t.shape):
    backstrides.add(i*(j-1))

  var iter_pos = t.offset

  ## Iterator loop
  for i in 0 ..< t.shape.product:

    ## Templating the return value
    when strider == IterKind.Values: yield t.data[iter_pos]
    elif strider == IterKind.Coord_Values: yield (coord, t.data[iter_pos])
    elif strider == IterKind.MemOffset: yield iter_pos
    elif strider == IterKind.MemOffset_Values: yield (iter_pos, t.data[iter_pos])

    ## Computing the next position
    for k in countdown(t.rank - 1,0):
      if coord[k] < t.shape[k]-1:
        coord[k] += 1
        iter_pos += t.strides[k]
        break
      else:
        coord[k] = 0
        iter_pos -= backstrides[k]
```
## in production: with init and next iterators templates

https://github.com/mratsim/Arraymancer/blob/3abe13b3a23d86536abe1151d4caba3e53f3dbe5/src/tensor/private/p_accessors.nim#L98-L131