xtensor-stack / xtensor

C++ tensors with broadcasting and lazy computing
BSD 3-Clause "New" or "Revised" License
3.33k stars 398 forks source link

question on view semantics #1350

Open cmauney opened 5 years ago

cmauney commented 5 years ago

Eigen recently introduced the seq and seqN functions such that:

using namespace Eigen;
Array<double, 10, 10> A, B, C
auto _i = seqN(1, fix<8>); // [1, 9)
auto _j = lastN(fix<4>); // [6, 10)
A(_i, _j) = 2.0 * B(_i, _j) + C(_i, _j);

with a simple operator overload, we can do something like (some kind of range checking implied)

// slightly simplified
template <typename Seq>
constexpr auto inline operator + (Seq s, int i) { return seqN(s.start() + i, s.size()); }
// ...
A(_i + 1, _j - 1 ) = B( _i + 1, _j ) + C(_i - 1, _j - 3); // _i + 1 == seqN(2, fix<8>)

We can do something similar in xtensor, with a caveat:

//overload as before, but with xt::range()
auto _i = range(1,9); auto _j = range(6,10);
view(A, _i + 1, _j - 1) = view(B, _i + 1, _j - 1) + view(C, _i - 1, _j - 3);

2 questions: -Is it possible to eliminate the prefix xt::view in expressions like this? If not, is it something that would be considered? -is this 'efficient' in xtensor? i.e. for the above example, are the ranges determined at compile-time? is the expression engine smart enough to end up with a single nested for loop, using offset indices?

Thanks!

JohanMabille commented 5 years ago

-Is it possible to eliminate the prefix xt::view in expressions like this? If not, is it something that would be considered?

It's not possible for now, views are built through calls to free functions. There are three reasons for that:

Before I expose the third reason, notice that these are choices that we made an I'm not opposed to a change (one could argue that the tight coupling is worth of it regarding the nice syntax it brings, and we can choose to always return xview since it's the best compromise features / performances).

Now the third reason is about consistency. Having operator() able to return views does not remove the need of a view free function. Indeed the view free function accepts fewer arguments than the number of dimensions. In that case, it appends xt::xall objects to the list of slices until it reaches the number of dimensions. For instance:

xt::xarray<double> arr1
  {{1.0, 2.0, 3.0},
   {2.0, 5.0, 7.0},
   {2.0, 5.0, 7.0}};

std::cout << xt::view(arr1, 1) << std::endl;
// => { 2., 5., 7. }
// Equivalent to 
std::cout << xt::view(arr1, 1, xt::all()) << std::endl;

On the other hand, the access operator() also accepts fewer arguments than the number of dimensions. In that case, it prepends 0 to the list of arguments until it reaches the number of dimension:

std::cout << arr1(2) << std::endl;
// Equivalent to std::cout << arr1(0, 2) << std::endl;

This behavior cannot be changed because it gives the guarantee that for any expressions with compatible shapes (w.r.t. the broadcasting rules ), the following is true:

(a + b)(i0, ...., iN) = a(i0, ..., iN) + b(i0, ...., iN);

So if you want to take a view on an expression by specifying integer slices only, you need the view free function. I would find it confusing to be able to build a view from operator() for some use cases, and need to call the free function view for other use cases.

SylvainCorlay commented 5 years ago

Yeah, also, making the choice statically depend on the number of argument can only work if dimensionality is a compile-time constant. Unlike eigen, xtensor support dynamically dimensioned arrays.

JohanMabille commented 5 years ago

-is this 'efficient' in xtensor? i.e. for the above example, are the ranges determined at compile-time? is the expression engine smart enough to end up with a single nested for loop, using offset indices?

As a complement to Sylvain's answer, even if the bounds of the range are not handled at compile time, the assignment of a view is optimized enough to produce a singled nested loop, where the assignment of contiguous parts of the expression can be accelerated with simd instructions.