Fast 4D argmax

mkitti · April 5, 2021, 4:53pm

How about using

julia_assign("A", A)
julia_eval("whichmax5Dt(A)", "Julia")

This should help us evaluate transfer issues.

jroon · April 5, 2021, 5:26pm

@mkitti - great suggestion - here is my R code and results now:

# Make an array
n = 20 ; nx = 400 ; v = 40 ; d = 3 
array4d <- array( runif(n*nx*v*d ,-1,0),
                  dim = c(n, nx, v, d) )
#Assign to Julia
julia_assign("A", array4d)

microbenchmark:: microbenchmark(
  res1 <- which_max_4D(array4d),
  res2 <- julia_eval("whichmax4D1(A)", "R"),
  res3 <- julia_eval("whichmax4D2(A)", "R"),
  res4 <- julia_eval("whichmax4Dmap(A)", "R"),
  res5 <- julia_eval("whichmax4Dreduce(A)", "R"),
  res6 <- julia_eval("whichmax4Dtullio(A)", "R"),
  res7 <- julia_eval("whichmax4D5_leandro(A)", "R"),
  res8 <- julia_eval("whichmax4D5_elrod(A)", "R"),
  res9 <- julia_eval("whichmax4D6(A)", "R"),
  res10 <- whichmaxC(array4d), times=500)

Wow!!

I also ran a benchmark on the julia_assign command:

Looks like you are right @mkitti - the slow part is moving the data from R to Julia

Edit: I ran it on an even bigger array also:

Non-Contradiction · April 5, 2021, 6:44pm

The performance here is indeed dominated by array transformation. Typically, the time can be decomposed into three parts.

JuliaCall interface function time cost (julia_call on the order of 10 microsends, julia_command on the order of 100 microseconds).
data transfer time cost, which will be on the order of milliseconds for a large data set like this.
julia function time cost.

stillyslalom · April 5, 2021, 7:13pm

It should be possible to wrap the raw bits of the R array instead of converting/copying - this is already done here in RCall.jl for R → Julia arrays.

Non-Contradiction · April 5, 2021, 8:14pm

Thanks for the information.
The current implementation of JuliaCall is on the “safe” side, so users have little to worry about. I have thought about having something like julia_unsafe_call which does not do things like automatic type conversion/data copying. I will think more about the idea when I have more time on this.
Another thing is that JuliaCall currently uses RCall.rcopy for data copying. I just noticed that the performance here seems to be magnitude slower than copy in julia. And in this case, the rcopy logic should be here: RCall.jl/base.jl at 0b5e9f119008d3a8471833566745311df79101ff · JuliaInterop/RCall.jl · GitHub
I’m wondering whether the performance for rcopy in this case can also be improved.

jroon · April 5, 2021, 8:34pm

That all sounds awesome! To give you my perspective on how I’d like to use JuliaCall → I’m fluent enough in R to be effective at getting my work tasks done quickly and dabble with making packages etc (nothing released yet), I expect it would take me a while to get to that level in Julia. But I would like to splash in snippets of Julia here and there were R is particularly slow, as opposed to writing large functions outright. I know many people use Rcpp in this way, but I’ve always found C++ to be utterly impenetrable for me. Learning Julia seems more achievable for me! So from that point of view, achieving blistering speed via a julia_unsafe_call for short functions sounds very appealing, but, I do understand why its something you need to think carefully about also! In any case JuliaCall is already an awesome package - thank you!

I must say, the response to my query here has been great and I’m blown away at all the options presented, and speeds that can be achieved - thank you all for the time you have put into this!

mkitti · April 6, 2021, 1:06am

To elaborate on the unsafe aspects, see the comment here:

github.com

JuliaInterop/RCall.jl/blob/0b5e9f119008d3a8471833566745311df79101ff/src/methods.jl#L43-L57


      
          """
          Represent the contents of a VectorSxp type as a `Vector`.
          
          
This does __not__ copy the contents.  If the argument is not named (in R) or
          otherwise protected from R's garbage collection (e.g. by keeping the
          containing RObject in scope) the contents of this vector can be modified or
          could cause a memory error when accessed.
          
          
The contents are as stored in R.  Missing values (NA's) are represented
          in R by sentinels.  Missing data values in RealSxp and CplxSxp show
          up as `NaN` and `NaN + NaNim`, respectively.  Missing data in IntSxp show up
          as `-2147483648`, the minimum 32-bit integer value.  Internally a `LglSxp` is
          represented as `Vector{Int32}`.  The convention is that `0` is `false`,
          `-2147483648` is `NA` and all other values represent `true`.
          """

The asymmetry of the transfers is immediately striking to me. I may be mistaken, but we do not see a lot of overhead of transferring the result from Julia to R. This suggests a strategy of starting with the variable on the Julia side.

Perhaps @randy3k or @simonbyrne could comment from the Rcall.jl perspective. I suspect julia_[unsafe]_call might not be that unsafe if julia_call could protect the data from R’s garbage collection. Some combination of robject, RCall.protect, and unsafe_array might work.