Search code examples
rmatlabbitwise-operatorshilbert-curve

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


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.


Solution

  • 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.

    enter image description here

    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')
    

    enter image description here

    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