How to translate bitwise functions from matlab/C to R? Particular case: Hilbert curve algorithm

Question

I am trying to translate a script written in matlab to R. The script maps 1D coordinates to 2D coordinates based on the Hilbert curve.

There is a line in the script which I am not sure how to translate into R:

ry = mod ( bitxor ( uint8 ( t ), uint8 ( rx ) ), 2 )

I think there is an R package with the bitxor() function, but not sure what to do about uint8().

Help appreciated!

The full matlab script can be found here:

https://people.sc.fsu.edu/~jburkardt/m_src/hilbert_curve/d2xy.m

The rot() function called in the script is here:

https://people.sc.fsu.edu/~jburkardt/m_src/hilbert_curve/rot.m

C versions can be found here:

https://en.m.wikipedia.org/wiki/Hilbert_curve

Some background, in case it's of interest: I am an amateur coder. Normally I only write programs where I understand the flow of logic from one line of code to the next. In this case I do not understand the logic, but I know what I want it to do and want it badly enough to go ahead rather blindly with this task.

In particular I have no clue what the bitxor() and unint8() functions are doing, although I understand what an xor logic gate is in principle.

I won't complain if some kind soul out there translates the whole script.

Knowing what _actual_ problem you're trying to solve would be much better but since you haven't told us that: https://github.com/hrbrmstr/hilbert — hrbrmstr, Jan 15 '17 at 03:29
@hrbrmstr The actual problem is to visualize a time series as a picture. For example, see what red, blue, brown, white noise looks like in 2D. Maybe a song or speech would be interesting to look at. I think your code will work for this, but I am having trouble downloading it. When I follow your instructions I get an error message telling me to install Rtools 3.4 from http://cran.r-project.org/bin/windows/Rtools/. But I cannot figure out how to do this. — ben, Jan 15 '17 at 04:03
I suggest do not use the program for practical purposes since it is very inefficient if you search you can find more efficient vectorized programs such as http://codegolf.stackexchange.com/a/102094/62451 — rahnema1, Jan 15 '17 at 04:52

Sathish · Accepted Answer · 2017-01-15T10:10:00.877

Matlab to R

# start d2xy    
d2xy <- function (m, d)
{
  m <- as.integer(m)
  d <- as.integer(d)

  n <- 2^m
  x <- 0
  y <- 0
  t <- d
  s <- 1

  while ( s < n ){
    rx <- floor ( t / 2 ) %% 2 
    if ( rx == 0 ){
      ry <- t %% 2
    } else {
      ry <- bitwXor(as.integer(t), as.integer(rx)) %%  2
    }

    xy <- rot ( s, x, y, rx, ry )
    x <- xy['x'] + s * rx
    y <- xy['y'] + s * ry
    t <- floor ( t / 4 )
    s <- s * 2
  }

  return(c(x = x, y = y))      
}
# end d2xy

# start rot
rot <- function(n, x, y, rx, ry) 
{
  n <- as.integer(n)
  x <- as.integer(x)
  y <- as.integer(y)
  rx <- as.integer(rx)
  ry <- as.integer(ry)

  if ( ry == 0 ){
    if ( rx == 1 ){
      x <- n - 1 - x
      y <- n - 1 - y
    }
    t <- x
    x <- y
    y <- t
  }

  return(c(x = x, y = y))
}
# end rot

Testing above functions in R

# vectorize our translated R function
d2xy_R <- Vectorize(d2xy, c('m', 'd'))    
rm(d2xy)

comparing matlab to R translated code with matlab functions

set.seed(1L)
m <- 2
d <- 5
xx <- runif(n = m*d, min = 0, max = 1)
mat_R <- d2xy_R(m = m, d = 1:d)
mat_R
#     [,1] [,2] [,3] [,4] [,5]
# x.x    1    1    0    0    0
# y.y    0    1    1    2    3

compare the mat_R output with matlab output. Both are same, hence no problem in translation.

mat_R <- mat_R + 1
coord2D_R <- matrix(xx[mat_R], nrow = m, ncol = d)
rownames(coord2D_R) <- c('x', 'y')

coord2D_R
#        [,1]      [,2]      [,3]      [,4]      [,5]
# x 0.3721239 0.3721239 0.2655087 0.2655087 0.2655087
# y 0.2655087 0.3721239 0.3721239 0.5728534 0.9082078

plot hilbert curve

set.seed(1L)
m <- 2
d <- 50
xx <- runif(n = m*d, min = 0, max = 1)
mat_R <- d2xy_R(m = m, d = 1:d)
mat_R <- mat_R + 1
coord2D_R <- matrix(xx[mat_R], nrow = m, ncol = d)
rownames(coord2D_R) <- c('x', 'y')
plot(t(coord2D_R), type = 'l', col = 'red')

Compare matlab and R translated codes with @hrbrmstr's github hilbert package

get hilbert.cpp file from hrbrmstr github hilbert package

library('Rcpp')
sourceCpp("hilbert.cpp") # compile C++ functions in hilbert.cpp file
d2xy_Rcpp <- d2xy
rm(d2xy)

mat_Rcpp <- matrix(nrow = m, ncol = d)
rownames(mat_Rcpp) <- c('x', 'y')

for(i in seq_len(d)){   # for loop is introduced, because unlike the R translated code, the Rcpp function is not vectorized
  xy <- d2xy_Rcpp(n = m, d = i)
  mat_Rcpp['x', i] <- xy['x']
  mat_Rcpp['y', i] <- xy['y']
}

mat_Rcpp
#      [,1] [,2] [,3] [,4] [,5]
# [1,]    0    1    1    0    0
# [2,]    1    1    0    0    1

compare mat_Rcpp output with mat_R and matlab outputs. It does not match with them, so there may be bugs in this package or the provided matlab code has issues with it.

mat_Rcpp <- mat_Rcpp + 1
coord2D_Rcpp <- matrix(xx[mat_Rcpp], nrow = m, ncol = d)
rownames(coord2D_Rcpp) <- c('x', 'y')

coord2D_Rcpp
#        [,1]      [,2]      [,3]      [,4]      [,5]
# x 0.2655087 0.3721239 0.3721239 0.2655087 0.2655087
# y 0.3721239 0.3721239 0.2655087 0.2655087 0.3721239

Benchmark matlab to R translated code with hrbrmstr's hilbert package

library('microbenchmark')
set.seed(1L)
m <- 2
d <- 5
xx <- runif(n = m*d, min = 0, max = 1)
microbenchmark(d2xy_R(m = m, d = d),      # matlab to R translation
               d2xy_Rcpp(n = m, d = d),   # @hrbrmstr - hilbert github package
               times = 100000)

# Unit: microseconds
#                    expr     min      lq       mean  median      uq       max neval
# d2xy_R(m = m, d = d)    169.382 177.534 192.422166 180.252 184.780 94995.239 1e+05
# d2xy_Rcpp(n = m, d = d)   2.718   4.530   7.309071   8.606   9.512  2099.603 1e+05

(a) why a list for the return? (b) it'd be interesting to compare that to github.com/hrbrmstr/hilbert — hrbrmstr, Jan 15 '17 at 03:48
forgot to add that it's a great translation example for operations like this. shld be a nice ref for others who are looking to do similar things with other algorithms. I'd prbly also add an `L` to the end of all the bare numerics to ensure they keep as integers for the speediest of maths ops. — hrbrmstr, Jan 15 '17 at 03:55
@Sathish Yes of course I am very grateful. One thing: After your edits, especially the sanity check, I am getting error:all(sapply(as.list(match.call()), is.integer)) is not TRUE. Here is what I am doing: ord<-2 dim2D<-2^ord nCells<-dim2D^2 xx<-runif(n = nCells, min = 0, max = 1) coord2D_list<-list() mat<-matrix(nrow=dim2D,ncol=dim2D) for(i in 1:nCells) { coord2D<-d2xy(ord,i) print(i) print(coord2D) mat[coord2D$x,coord2D$y]<-xx[i] } — ben, Jan 15 '17 at 04:48
@Sathish I have set both i and ord to integer (using as.integer()), and I included the set.seed(1L) at the top, but still get error. Another question: Note that the matrix is bordered with NA. Shouldn't d2xy() output coordinates for all cells of a 2^ord-by-2^ord matrix? — ben, Jan 15 '17 at 05:10
@hrbrmstr I think you may have a bug in your hilbert package. Please see my edited answer and please confirm. Thanks — Sathish, Jan 15 '17 at 08:34
cld very well be a bug but it's not my code. it's the copy/pasted C code from a random link from the OP. — hrbrmstr, Jan 15 '17 at 09:39

How to translate bitwise functions from matlab/C to R? Particular case: Hilbert curve algorithm

1 Answers1