xtensor-stack / xtensor-r

R bindings for xtensor
BSD 3-Clause "New" or "Revised" License
86 stars 15 forks source link

Is eye() specific to row major layouts? #90

Open DavisVaughan opened 5 years ago

DavisVaughan commented 5 years ago

I can't tell if this is expected behavior of eye() or not. Is there a column major option for it?

I dont really see anything here about row/column major, so I honestly have no idea https://github.com/QuantStack/xtensor/blob/master/include/xtensor/xbuilder.hpp#L305

// [[Rcpp::depends(xtensorrr)]]
// [[Rcpp::plugins(cpp14)]]

#include <xtensor-r/rarray.hpp>
#include <Rcpp.h>
using namespace Rcpp;

// [[Rcpp::export]]
SEXP rray_eye_cpp(const std::vector<std::size_t> shape, int k = 0) {
  xt::rarray<int> res = xt::eye<int>(shape, k);
  return res;
}
Rcpp::sourceCpp("~/Desktop/test.cpp")

# looks good
rray_eye_cpp(c(2, 2))
#>      [,1] [,2]
#> [1,]    1    0
#> [2,]    0    1

rray_eye_cpp(c(2, 2, 2))
#> , , 1
#> 
#>      [,1] [,2]
#> [1,]    1    0
#> [2,]    1    0
#> 
#> , , 2
#> 
#>      [,1] [,2]
#> [1,]    0    1
#> [2,]    0    1

# I think I would have expected
array(c(1,0,0,1,1,0,0,1), c(2,2,2))
#> , , 1
#> 
#>      [,1] [,2]
#> [1,]    1    0
#> [2,]    0    1
#> 
#> , , 2
#> 
#>      [,1] [,2]
#> [1,]    1    0
#> [2,]    0    1

Created on 2019-01-09 by the reprex package (v0.2.1.9000)

wolfv commented 5 years ago

that's definitely a bug it should behave the same way for row- vs col-major (your expectation seems correct)

wolfv commented 5 years ago

looking into this right now, so far i couldn't reproduce with "pure" xtensor. checking xtensor-r now!

wolfv commented 5 years ago

Ok, I looked into this and I think this is related to R's way of printing.

E.g. check out:

a[1,1,1] == 1
a[1,1,2] == 0
a[1,2,1] == 0
a[1,2,2] == 1
a[2,1,1] == 1
... all 0 ... 
a[2,2,2] == 1

and this is true and works. However, the way R prints multidimensional arrays is very awkward, and leads to this confusion.

What do you think @DavisVaughan ?

Sorry, again, that it took so long to get back to this.

wolfv commented 5 years ago

Ok, just thought about it for two more seconds, and i guess what you want is to have an eye function that is 1 where the first two dimensions are equal, so that

a[1,1,1] == a[2,2,1] == a[1,1,2] == a[2,2,2] == 1

Which I guess makes sense if one would do column major style broadcasting for things like matrix product etc, or anything else.

Now this is again leading to more fundamental questions ... wether we could / should support col-major "broadcasting". I fear right now we don't have the manpower for that and it's also not clear that it would bring a lot of advantages. On the other hand, given all the code we already have for row-major broadcasting, it should not be so far out to do it, but it adds some complexity to the overall code base...

DavisVaughan commented 5 years ago

Ah this looks heavily related to the prepend / append of new dimension conversation in #57.

Where you have a[1,2,2]=1 (Python and xtensor convention), in the R world we would have subset this as a[2,2,1]=1

To me, both mean "second element of first dimension, second element of second dimension, first element of third dimension"

I don't think it has to do with the way R prints. To be fair, because of these prepend/append of dimension differences (which bleeds into how things print), I think numpy array printing is highly confusing, with R printing being more intuitive.

About to get on a plane! Can chat more later today

DavisVaughan commented 5 years ago

The new equality you show above looks correct to me!

DavisVaughan commented 5 years ago

I think eye() needs an order/layout template parameter as well, as numpy.eye() has an order = 'C' argument. I'm fairly certain that's the fix here. https://docs.scipy.org/doc/numpy/reference/generated/numpy.eye.html