xtensor-stack / xtensor-r

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

Broadcasting to same dimension reorders #103

Closed DavisVaughan closed 5 years ago

DavisVaughan commented 5 years ago

I'm guessing this is just a column major thing, but I hit a weird case where broadcasting to the exact dimensionality that you are already at reorders the elements. Interestingly, broadcasting with an actual shape change works fine.

library(rray)

x <- rray(1:4, c(1, 2, 2, 1))
x
#> <vctrs_rray<integer>[,2,2,1][4]>
#> , , 1, 1
#> 
#>      [,1] [,2]
#> [1,]    1    2
#> 
#> , , 2, 1
#> 
#>      [,1] [,2]
#> [1,]    3    4

# Hmm - broadcast to same dimensions
rray_broadcast(x, c(1, 2, 2, 1))
#> <vctrs_rray<integer>[,2,2,1][4]>
#> , , 1, 1
#> 
#>      [,1] [,2]
#> [1,]    1    3
#> 
#> , , 2, 1
#> 
#>      [,1] [,2]
#> [1,]    2    4

# Behavior not seen when we broadcast to 
# other dimensions
rray_broadcast(x, c(1, 2, 2, 2))
#> <vctrs_rray<integer>[,2,2,2][8]>
#> , , 1, 1
#> 
#>      [,1] [,2]
#> [1,]    1    2
#> 
#> , , 2, 1
#> 
#>      [,1] [,2]
#> [1,]    3    4
#> 
#> , , 1, 2
#> 
#>      [,1] [,2]
#> [1,]    1    2
#> 
#> , , 2, 2
#> 
#>      [,1] [,2]
#> [1,]    3    4

Created on 2019-03-28 by the reprex package (v0.2.1.9000)

DavisVaughan commented 5 years ago

Here is a way to explore it without rray.

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

#include <iostream>
#include <xtensor-r/rarray.hpp>
#include "xtensor/xarray.hpp"
#include "xtensor/xio.hpp"
#include "xtensor/xview.hpp"
#include <Rcpp.h>
using namespace Rcpp;

// [[Rcpp::export]]
SEXP broadcast_test(SEXP x, SEXP arg) {
  xt::rarray<double> x_rarray(x);

  std::cout << "x_rarray:\n" << x_rarray << std::endl;

  std::vector<std::size_t> dim = Rcpp::as<std::vector<std::size_t>>(arg);

  xt::rarray<double> out = xt::broadcast(x_rarray, dim);

  std::cout << "out:\n" << out << std::endl;

  return out;
}
Rcpp::sourceCpp("~/Desktop/test.cpp")

x <- array(c(1, 2, 3, 4), c(1, 2, 2, 1))

x
#> , , 1, 1
#> 
#>      [,1] [,2]
#> [1,]    1    2
#> 
#> , , 2, 1
#> 
#>      [,1] [,2]
#> [1,]    3    4

broadcast_test(x, c(1, 2, 2, 1))
#> , , 1, 1
#> 
#>      [,1] [,2]
#> [1,]    1    3
#> 
#> , , 2, 1
#> 
#>      [,1] [,2]
#> [1,]    2    4

broadcast_test(x, c(1, 2, 2, 2))
#> , , 1, 1
#> 
#>      [,1] [,2]
#> [1,]    1    2
#> 
#> , , 2, 1
#> 
#>      [,1] [,2]
#> [1,]    3    4
#> 
#> , , 1, 2
#> 
#>      [,1] [,2]
#> [1,]    1    2
#> 
#> , , 2, 2
#> 
#>      [,1] [,2]
#> [1,]    3    4

# here is what the std output would look like
broadcast_test(x, c(1, 2, 2, 1))
x_rarray:
{{{{ 1.},
   { 3.}},
  {{ 2.},
   { 4.}}}}
out:
{{{{ 1.},
   { 2.}},
  {{ 3.},
   { 4.}}}}

broadcast_test(x, c(1, 2, 2, 2))
x_rarray:
{{{{ 1.},
   { 3.}},
  {{ 2.},
   { 4.}}}}
out:
{{{{ 1.,  1.},
   { 3.,  3.}},
  {{ 2.,  2.},
   { 4.,  4.}}}}
DavisVaughan commented 5 years ago

This problem is actually super obnoxious as I start to build other functions on top of rray_broadcast(). No pressure if you are super busy at the moment, but if you work on xtensor-r at any point I'd love to see this get high priority. I've been trying to debug, but can't seem to get far.

SylvainCorlay commented 5 years ago

Hey @DavisVaughan, thanks for the report!

There appears to be a legitimate issue in column-major mode for xt::broadcast for column-major layouts:

Screenshot from 2019-03-30 01-54-24

Will be investiating.

DavisVaughan commented 5 years ago

Thanks!

SylvainCorlay commented 5 years ago

xt::broadcast expression access appears ok, but iteration and assignment look very wrong (in column major order).

SylvainCorlay commented 5 years ago

it seems that the issue is with the linear assigner.