For discussion: use dedicated shape objects

[ukoethe]

This PR offers an implementation that replaces the current std::array and std::vector shapes with a dedicated shape class tiny_array. The PR is mainly intended for discussion, although it might eventually be merged if people happen to like the proposed shape class 😃 . I successfully executed the full xtensor test suite using Visual Studio, but did not yet test much on other compilers. Changes come in three groups:

Preliminaries

The PR adds headers xconcepts.hpp (concept checking and type inference), xmathutil.hpp (additional mathematical functions) and xtags.hpp (simple keyword argument support). Concept checking can now be done via

template <class T, 
          XTENSOR_REQUIRE<std::is_arithmetic<T>::value>>
T foo(T t);

which is arguably more readable than enable_if_t applied to the return type.

Shape Class

Files xtiny.hpp and test_xtiny.cpp implement and test the classes tiny_array (arrays that own their memory) and tiny_array_view (arrays that don't own their memory) along with a large set of arithmetic, logical, and algebraic functions. tiny_array<size_t, 3> corresponds to std::array<size_t, 3>, and tiny_array<size_t, runtime_size> corresponds to std::vector<size_t>. The shape classes are designed to resemble the behavior of built-in arithmetic types. In particular, operations usually return new objects, rather than modifying existing ones in-place. This is key to the unification of static and dynamic shapes. IMHO, this design decision won't cause performance problems, because shape calculations are rarely done in the inner loop.

It is open for debate if it would be prefarable to implement mathematical functions for tiny_array in terms of the existing xexpression machinery. A possible advantage of the present design is that the compiler will unroll loops automatically for static shapes. (As an aside, I had to add XTENSOR_REQUIRE<xexpression_concept<E>::value> to some functions in xmath.hpp because these templates matched too greedily and thus dominated the corresponding functions in xtiny.hpp.)

Changes to the Existing Code Base

In the original code, declarations of shape types are spread all over the place (I found >30 occurences in the headers and another ~60 in the tests). I first implemented a single point of reference for these declarations by defining aliases stat_shape and dyn_shape in xutils.hpp and replacing hard-coded declarations of std::array and std::vector as appropriate. I then added overloads of the required utility functions and metafunctions for tiny_array and finally replaced the types in stat_shape and dyn_shape with tiny_array. Enjoy!

A major advantage of the present design is that xarray and xtensor essentially become the same class -- all differences can potentially be handled by tiny_array.

[wolfv]

Argh, I can't comment on the lines above, but there is still a std::enable_if_t<!detail::is_array<SC>::value... which obviously should be replaced. Maybe we can have something like detail::is_static_container or detail::is_static_size?
This currently prevents building of the xadapt tests with GCC.

[wolfv]

Sorry, I just discovered your new overload further down. Unfortunately it doesn't seem to be picked up. I'd love to investigate more, but lack the time (it's on Ubuntu 16.04 / GCC 5.4).

[wolfv]

I think in terms of naming i would prefer xshape and then decide on it being static vs. dynamic based on the template arguments. Is it necessary to have dyn_shape vs stat_shape?

[wolfv]

Ie. xshape<3> == stat_shape, and xshape() == dyn_shape...

[wolfv]

On that note, I also think we shouldn't use the size_t template everywhere. It would be much better to define the xshape<xt::index_t> somewhere and then have it consistent across the library and only one place to change it.

[ukoethe]

I agree that xshape<3> and xshape<> are better alias types in the long run. During the experimentation phase, the present definition has the advantage of making the old and new shape types exchangeable.

[wolfv]

I know it's nitpicky, but it would be cool to have a similar style everywhere. We don't use these indented triple-slash comments. :)

[ukoethe]

Triple-slashs are simply doxygen's C++ variant of /** for one-line documentation. I don't see the necessity to enforce a convention here, but would change the code if you insist.

[ukoethe]

Moreover, I consider the code more readable when the documentation is indented more than the code.

[wolfv]

And we always put newlines in front of brackets

[ukoethe]

So do I, with this single exception: I usually don't indent at namespace level. If you prefer, I can indent the code.

[wolfv]

I think it's a really great addition to xtensor! Thanks for contributing.

However, I wonder wether it would be better to remove some features for now, and add them at a later stage?

If we concentrate on the TinyArray being a shape type, do we really need all the math functions defined for it? E.g. the trigonometric stuff etc.

Also, it might be confusing to have another kind of view in xtensor which doesn't adhere to the same interface as the existing views.
In general I guess it's easier to add features than to remove them ...

For example, I am not sure if we really need all_greater and related functions ... with other xexpressions you'd generally do xt::all(a > b) to get the same effect. Johan and I were quickly discussing the idea to make the TinyArray into an xexpression itself by inheriting from xcontainer ... I think conceptually it would be cool, and we'd save all the declarations for math etc. but maybe it also adds some overhead.

What do you think?

I hope this doesn't come across as too negative. I think this is really important work and personally I really want to see it in xtensor.

[DerThorsten]

Is it possible to make the shape an expression?
The shape itself is somehow needed within the expression system right?
I guess this would give an infinite template recursion.
How would you bootstrap that?
Seems like a nice meta-programming challenge =).

[wolfv]

Luckily we can define the shape type of the shape type to be something else ;)

So xshape-container shape could be a std::array<xt::index_t, 1> container! (Which is easy to assume as shapes will always be 1D). The question is if this carries any overhead (which I can't answer).

[JohanMabille]

Each xexpression has a shape type. In order to make the shape an xexpression, we only need to give it a particular shape type that would not be an xexpression (that could be an std::array with a single element for instance).

Edit: @wolfv you beat me to it!

[DerThorsten]

Ok, but if the ShapeType is an expression, and the shape of that expression is again something very special, then whats the point of having the ShapeType shape.
So I guess shapes would need a different expression mechanism (aka a shape CRTP base).
Or one just implements all the math/expressions for the shapes by hand.

[JohanMabille]

You need to perform mathematical operations on the shapes, not on the size of the shape (which would become the shape of the shape). Making the shape_type an xexpression is just a way to avoid code duplication and reuse the implementation of the math functions. The special type used as shape_type of the shape has no importance, it's only meant to stop the template recursion.

[ukoethe]

If we concentrate on the TinyArray being a shape type, do we really need all the math functions defined for it? E.g. the trigonometric stuff etc.

I was of the same opinion initially (a decade or so ago 😄), but functions make all the difference. The most important lesson I learned over the years is "Users love convenience functions." Just consider

dot(point, strides)  // memory address offset
prod(shape)          // number of array elements
norm(shape)          // length of the diagonal

Since I always forgot which functions were available and which weren't, I adopted the easy-to-remember rule "support everything from <cmath> plus some numpy-inspired array-specific stuff".

Also, it might be confusing to have another kind of view in xtensor which doesn't adhere to the same interface as the existing views.

True, but tiny_array_view is needed for the subarray function. For example

auto backstrides = stride_type::unit_vector(0);
backstrides.subarray<1, 3>() = cumprod(shape.subarray<0,2>());

(ok, this doesn't special-case singleton dimensions ATM).

For example, I am not sure if we really need all_greater and related functions ... with other xexpressions you'd generally do xt::all(a > b) to get the same effect.

all_greater and xt::all(a > b) are very different here: in order to conform to the std::vector-API, I implemented a > b as a lexicographic comparison returning bool. I also considered implementing it element-wise, but came to the conclusion that lexicographic was more useful (e.g., one can put tiny_arrays into a std::map) and more consistent (would any user expect a == b to be element-wise for a shape class?). OTOH, if you have good arguments, the design is easy to change.

Johan and I were quickly discussing the idea to make the TinyArray into an xexpression itself by inheriting from xcontainer.

I don't know the details of xexpression sufficiently to judge if this would be a good idea. It might avoid some code duplication, but would make the design harder to understand. In the particular case of a shape class, code duplication is not a very big concern because the arithmetic and algebraic functions are unlikely to ever change once implemented. I also had some trouble with circular dependencies when I tried to incorporate the contents of xmathutil.hpp into xmath.hpp (xtiny.hpp than included xmath.hpp which includes xutils.hpp which in turn was supposed to include xtiny.hpp). These problems are certainly solvable, but I'm not sure if the advantages are worthwhile.

[wolfv]

I am definitly not arguing against having all the syntactic sugar, I think it's awesome. But I think it will be much easier to merge a bunch of small commits. And I think we need to be careful to integrate the ideas from xtiny into the overall framework. It might be tough to add a lot of functions now, only to decide to e.g. rename some of them later.

Couldn't we do a PR with only the tinyarray and the arithmetic functions first, and then do separate PR's for views, and more functionality? At least that's what I think would be easier to review and figure everything out.
But maybe it's also worth it to hold back a little, as Sylvain and Johan might have some other feedback.

[wolfv]

I'm running into some weird compiler bug it seems:

[ 28%] Building CXX object test/CMakeFiles/test_xtensor.dir/test_xio.cpp.o
In file included from /shape/xtensor/include/xtensor/xview.hpp:24:0,
                 from /shape/xtensor/include/xtensor/xio.hpp:22,
                 from /shape/xtensor/test/test_xio.cpp:18:
/shape/xtensor/include/xtensor/xview_utils.hpp: In instantiation of ‘struct xt::detail::view_temporary_type_impl<double, xt::tiny_array<long unsigned int, -1>, (xt::layout_type)1, int, xt::xall<long unsigned int> >’:
/shape/xtensor/include/xtensor/xview_utils.hpp:98:76:   required from ‘struct xt::view_temporary_type<xt::xarray_container<xt::uvector<double, std::allocator<double> >, (xt::layout_type)1, xt::tiny_array<long unsigned int, -1> >, int, xt::xall<long unsigned int> >’
/shape/xtensor/include/xtensor/xview_utils.hpp:102:79:   required by substitution of ‘template<class E, class ... SL> using view_temporary_type_t = typename xt::view_temporary_type::type [with E = std::decay<xt::xarray_container<xt::uvector<double, std::allocator<double> >, (xt::layout_type)1, xt::tiny_array<long unsigned int, -1> >&>::type; SL = {int, xt::xall<long unsigned int>}]’
/shape/xtensor/include/xtensor/xview.hpp:37:77:   required from ‘struct xt::xcontainer_inner_types<xt::xview<xt::xarray_container<xt::uvector<double, std::allocator<double> >, (xt::layout_type)1, xt::tiny_array<long unsigned int, -1> >&, int, xt::xall<long unsigned int> > >’
/shape/xtensor/include/xtensor/xsemantic.hpp:39:82:   required from ‘class xt::xsemantic_base<xt::xview<xt::xarray_container<xt::uvector<double, std::allocator<double> >, (xt::layout_type)1, xt::tiny_array<long unsigned int, -1> >&, int, xt::xall<long unsigned int> > >’
/shape/xtensor/include/xtensor/xsemantic.hpp:201:11:   required from ‘class xt::xview_semantic<xt::xview<xt::xarray_container<xt::uvector<double, std::allocator<double> >, (xt::layout_type)1, xt::tiny_array<long unsigned int, -1> >&, int, xt::xall<long unsigned int> > >’
/shape/xtensor/include/xtensor/xview.hpp:74:11:   required from ‘class xt::xview<xt::xarray_container<xt::uvector<double, std::allocator<double> >, (xt::layout_type)1, xt::tiny_array<long unsigned int, -1> >&, int, xt::xall<long unsigned int> >’
/shape/xtensor/include/xtensor/xview.hpp:744:13:   required from ‘auto xt::detail::make_view_impl(E&&, std::index_sequence<I ...>, S&& ...) [with E = xt::xarray_container<xt::uvector<double, std::allocator<double> >, (xt::layout_type)1, xt::tiny_array<long unsigned int, -1> >&; long unsigned int ...I = {0ul, 1ul}; S = {int, xt::xall_tag}; std::index_sequence<I ...> = std::integer_sequence<long unsigned int, 0ul, 1ul>]’
/shape/xtensor/include/xtensor/xview.hpp:759:38:   required from ‘auto xt::view(E&&, S&& ...) [with E = xt::xarray_container<xt::uvector<double, std::allocator<double> >, (xt::layout_type)1, xt::tiny_array<long unsigned int, -1> >&; S = {int, xt::xall_tag}]’
/shape/xtensor/test/test_xio.cpp:59:40:   required from here
/shape/xtensor/include/xtensor/xview_utils.hpp:88:38: internal compiler error: in tsubst, at cp/pt.c:12186
             using type = xarray<T, L>;
                                      ^

When I change the T to some hard coded value (e.g. double) it seems to work.

Edit: this is on GCC 5.4 / Ubuntu 16.04. I'll try later with GCC 7 on Fedora 26.

[wolfv]

With clang 3.8 it builds if I remove the constexpr on line 1723 in xtiny.hpp as the return type is not a LiteralType.
I also commented out the traits tests in test_xtiny as they failed to link.

[wolfv]

GCC 7 seems to work with the exception of xadapt. I think I found the issue: The templates for static shapes / xtensor are using std::size_t N while the tiny_array uses an int N and therefore no overload is found. I can try to debug this tomorrow.

[ukoethe]

I wouldn't mind to get rid of tiny_array_view entirely (especially as this would also eliminate the need for tiny_array_base). However, I would not like to loose the subarray() function on the way. If someone comes up with a better way of distinguishing arrays that own their memory from those that don't, I'd be all ears. One possibility is to define

template <class T, bool owns_memory, int ... N>
class tiny_array;

but this doesn't look terribly elegant either. Edit: see #374 (comment) below for a better idea.

[ukoethe]

I can split the PR into several parts, but not before end of August (I'm heading off to vacation tomorrow 😀 ). And of course, once the tests succeed on all platforms, there are a ton of warnings waiting to be fixed. Moreover, if the PR gets accepted one might want to change some of the existing shape manipulation code to take advantage of the new syntactic sugar.

[ukoethe]

Another important point that needs to be discussed: right now, the shape's value_type is size_t. I'd strongly recommend to change this into a signed type such as std::ptrdiff_t, so that something like abs(shape1 - shape2) would work without unpleasant surprizes.

[ukoethe]

How about the following idea to eliminate tiny_array_view and tiny_array_base?

template <class T, bool owns_memory, int ... N>
class tiny_array_impl;

template <class T, int M=runtime_size, int ... N>
using tiny_array = tiny_array_impl<T, true, M, N...>;

template <class T, int M=runtime_size, int ... N>
using tiny_array_adaptor = tiny_array_impl<T, false, M, N...>;

template <class T, int M, int N>
using tiny_matrix = tiny_array_impl<T, true, M, N>;

[SylvainCorlay]

@ukoethe after reading your proposal, I have a better understanding of where you want to go.

I wanted to say that this is really sophisticated work. And it also shows that you guys really have a real mastery of the type of metaprogramming that is in play here.

I would like to propose a plan to move forward with these ideas. However, rather than integrating a massive change like this as a first PR, I would like to have a more progressive approach.

1) Static and dynamic shapes

One feature of both xtensor and in your work is that we both support dynamic and static dimensionality. There is a huge gain in performances from using statically-sized shapes and strides. The dynamic case is required for wrapping the containers from interpreted languages such as python (numpy) and julia, in xtensor-python and xtensor-julia.

A first way to reduce the gap between stack-allocated shapes and heap-allocated shapes would be to use a dynamically-sized container type that

• would have a stack-allocated buffer of constant size, used until a certain dimension N (e.g. N=4)
• a heap allocated buffer, that would only be used for higher dimensions.

Such a container would have the same performances as std::array for D<=N, (the most common case) and the same performances as std::vector for higher dimensions.

A good place to start the implementation of this new container would be to start from our uvector container. uvector differs from std::vector in that it does not initialize its values, and is used as a storage type in our container expressions (xarray, xtensor). Then, everywhere we use std::vector as a shape type by default, such as in xarray, we would use that new xvector. This change only would reduce the gap between the performances of xtensor and xarray.

2) Update to the promote_shape mechanism

Currently in xtensor, if we do auto e = a + b / d * sin(d) where all the members of the expressions are xexpressions, e is an un-evaluated expression. The shape type for e is

• std::array if all the members of the expression have shape type std::array
• std::vector if any of the members of the expression have another shape type.

Instead of having std::array as the only shape type considered "static", we would

• add a constexpr static boolean in our custom shape types specifying if it is static or dynamic.
• use a traits mechanism to handle the containers from the STL.

Then, downstream could use their custom shape types that would be recognized as static.

3) arithmetic and extra operations on shapes

• xtensor and xarray are just template aliases for specific shape types. We document how a use can create a custom expression (that would work with xtensor and xarray) with their own data and shapes containers... I think that it is a very powerful thing, and the only requirement here is that these types implement the standard interface for an STL container. Requiring more would limit the ability to define your own expressions...

• However, I get the need for arithmetic operations on shapes. A means to do this without duplicating all the math operations from xtensor would be for shape() to return an adaptor as an xexpression on the underlying shape type. That would remove a lot of code, since one would not need to re-implement all the arithmetic operations and features for the shape types, while they are already in xtensor.

4) Conclusion

Apologies for taking so much time to give you proper feedback. We should probably talk viva voce (over skype, hangout or appear.in) sometimes.

Let me know if you like this plan ! Items 1) and 2) and 3) can be done in independent PRs.

[DerThorsten]

to item 1)
We implemented such thing within OpenGM
https://github.com/DerThorsten/opengm/blob/master/include/opengm/datastructures/fast_sequence.hxx
But proper benchmarking needs to be done to ensure that it is as fast as std::array.
(The OpenGM impl. is C++ 98)

[ukoethe]

We should probably talk viva voce.

That would be good!

Rather than integrating a massive change like this as a first PR, I would like to have a more progressive approach.

I think it is beneficial for our discussion when you and others can look at the PR in its entirety. Actual integration - if agreed upon - can then proceed more gradually.

Static and dynamic shapes

The most important benefit of my design is that static and dynamic shapes have the same API and can be used interchangeably in generic code. In my experience, this leads to huge simplifications in functionality like promote_shape, dynamic ROI definition etc. Shape arithmetic in this PR is implemented such that static shapes are used whenever the resulting dimension can be determined at compile time, and dynamic shapes otherwise (this happens automatically, no additional template magic is needed).

heap allocated shapes

I'm not sure if allocation of the shape storage is really the bottleneck. According to my experience, the most critical operation is the transformation of index objects into memory addresses (i.e. the computation of the scalar product between index and strides), which regularly occurs in inner loops. The compiler can optimize this much much better if the dimension is known at compile time, but this optimization is independent of the allocation mechanism. If benchmarking experiments demonstrate that heap allocation is indeed too slow, I can implement the mixed allocation mode you proposed (stack for small N, heap for bigger N) or a custom allocator that calls new for batches of shape objects. Without such benchmarks, replacing heap allocation might be premature.

Then, downstream could use their custom shape types that would be recognized as static.

I'm proposing to replace std::vector and std::array in their role as shape types with tiny_array in the entire library, because tiny_array fully unifies static and dynamic shapes and provides a rich set of arithmetic and algebraic functions. I don't (yet) see the benefit of retaining std::vector and std::array and supporting additional custom shape types - this sounds unnecessarily complicated.

The dynamic case is (only) required for wrapping the containers from interpreted languages.

I found over the years that dynamic tensors are more generally useful. For example, functions which are separable over dimensions (e.g. Gaussian filters, distance transforms, data I/O) don't run faster if the dimensionality is known at compile time, whereas the benefits of dynamic dimensions are substantial: compilation times are greatly reduced and binaries get much smaller, because the compiler needs to instantiate these templates only once instead of repeatedly for every dimension.

A means to do this without duplicating all the math operations from xtensor...

In general, I agree that redundancy in code should be avoided. However, I'm not so sure about the present case:

• Using xexpressions for shape computation makes the inner workings of tiny_array a lot harder to understand. This may be a problem because many programmers don't master advanced template techniques, which means that they have a hard time interpreting error messages for even trivial mistakes and are less likely to contribute back to the code base. 1200 extra lines (minus the required adaptor code) of easy arithmetic might be a preferable trade-off.
• Redundancy is not that big a problem if the code is unlikely to ever change, as I expect to be true for shape arithmetic.
• I'm not sure if the same arithmetic rules apply or should apply to tensors and shapes (e.g. shape1 < shape2 implements lexicographic comparison returning bool, to be consistent with std::vector).

[SylvainCorlay]

We should probably talk viva voce.

That would be good!

I am coming back from JupyterCon on wednesday. Both I and Johan would be available pretty much anytime Thurday or Friday, or next week.

Rather than integrating a massive change like this as a first PR, I would like to have a more progressive approach.

I think it is beneficial for our discussion when you and others can look at the PR in its entirety. Actual integration - if agreed upon - can then proceed more gradually.

Awesome, I think that we are on the same page.

[...]

I'm not sure if allocation of the shape storage is really the bottleneck. According to my experience, the most critical operation is the transformation of index objects into memory addresses (i.e. the computation of the scalar product between index and strides), which regularly occurs in inner loops. The compiler can optimize this much much better if the dimension is known at compile time, but this optimization is independent of the allocation mechanism. If benchmarking experiments demonstrate that heap allocation is indeed too slow, I can implement the mixed allocation mode you proposed (stack for small N, heap for bigger N) or a custom allocator that calls new for batches of shape objects. Without such benchmarks, replacing heap allocation might be premature.

I agree with you about the cost of indexed-based indexing, but the iterators that we have in xtensor don't rely on such inner products. So in the case of iterator-based indexing, the heap allocation becomes dominant.

Then, downstream could use their custom shape types that would be recognized as static.

I'm proposing to replace std::vector and std::array in their role as shape types with tiny_array in the entire library, because tiny_array fully unifies static and dynamic shapes and provides a rich set of arithmetic and algebraic functions. I don't (yet) see the benefit of retaining std::vector and std::array and supporting additional custom shape types - this sounds unnecessarily complicated.

We have a very strong use-case for supporting user-provided shape types, for the adaptors from other data structures. (e.g. in xtensor-python, the inner shape type is a buffer adaptor). However, I hear your need for richer arithmetic on shapes, so I have a proposal which should address this usecase.

• only require the inner shape type to have the API of a standard container.
• have shape() return a in-place adaptor as an xexpression on those.

So that user-defined expressions adapting existing data structures don't need to recursively define expressions for the shapes etc...

The dynamic case is (only) required for wrapping the containers from interpreted languages.

I found over the years that dynamic tensors are more generally useful. For example, functions which are separable over dimensions (e.g. Gaussian filters, distance transforms, data I/O) don't run faster if the dimensionality is known at compile time, whereas the benefits of dynamic dimensions are substantial: compilation times are greatly reduced and binaries get much smaller, because the compiler needs to instantiate these templates only once instead of repeatedly for every dimension.

👍 I agree with that 100%

A means to do this without duplicating all the math operations from xtensor...

In general, I agree that redundancy in code should be avoided. However, I'm not so sure about the present case:

Using xexpressions for shape computation makes the inner workings of tiny_array a lot harder to understand. This may be a problem because many programmers don't master advanced template techniques, which means that they have a hard time interpreting error messages for even trivial mistakes and are less likely to contribute back to the code base. 1200 extra lines (minus the required adaptor code) of easy arithmetic might be a preferable trade-off.
Redundancy is not that big a problem if the code is unlikely to ever change, as I expect to be true for shape arithmetic.

I think that the adaptor approach proposed above bridges the gap by avoiding the duplication of the code, while proving a rich API for shape types and making shape type interoperable with expressions.

I'm not sure if the same arithmetic rules apply or should apply to tensors and shapes (e.g. shape1 < shape2 implements lexicographic comparison returning bool, to be consistent with std::vector).

Gotcha. I don't think that is a big issue.

[SylvainCorlay]

Closing this one which, which was split into multiple discussions.