如何将按位函数从 matlab/C 转换为 R？特例：希尔伯特曲线算法

Posted 2023-04-17

【中文标题】如何将按位函数从 matlab/C 转换为 R？特例：希尔伯特曲线算法

【英文标题】：How to translate bitwise functions from matlab/C to R? Particular case: Hilbert curve algorithm

【发布时间】：2017-05-30 04:02:54

【问题描述】：

我正在尝试将用 matlab 编写的脚本翻译成 R。该脚本根据希尔伯特曲线将 1D 坐标映射到 2D 坐标。

脚本中有一行我不知道如何翻译成 R：

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

我认为有一个带有 bitxor() 函数的 R 包，但不确定如何处理 uint8()。

帮助表示赞赏！

完整的matlab脚本可以在这里找到：

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

脚本中调用的rot()函数在这里：

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

C 版本可以在这里找到：

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

一些背景，以防感兴趣：我是一名业余编码员。通常我只编写我理解从一行代码到下一行代码的逻辑流程的程序。在这种情况下，我不明白其中的逻辑，但我知道我想要它做什么，并且非常希望它能够盲目地继续这项任务。

特别是我不知道 bitxor() 和 unint8() 函数在做什么，尽管我理解 xor 逻辑门的原理。

如果有好心人翻译整个剧本，我不会抱怨。

【问题讨论】：

了解您要解决的实际问题会好得多，但既然您还没有告诉我们：github.com/hrbrmstr/hilbert

@hrbrmstr 实际问题是将时间序列可视化为图片。例如，查看红色、蓝色、棕色、白色噪声在 2D 中的样子。也许一首歌或演讲会很有趣。我认为您的代码适用于此，但我无法下载它。当我按照您的说明进行操作时，我收到一条错误消息，告诉我从cran.r-project.org/bin/windows/Rtools 安装 Rtools 3.4。但我不知道该怎么做。

我建议不要将该程序用于实际目的，因为它非常低效，如果您搜索您可以找到更有效的矢量化程序，例如codegolf.stackexchange.com/a/102094/62451

【参考方案1】：

Matlab 转 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

在 R 中测试上述函数

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

使用 matlab 函数将 matlab 与 R 翻译代码进行比较

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

将mat_R 输出与matlab 输出进行比较。两者都是一样的，因此翻译没有问题。

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

绘制希尔伯特曲线

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

将 matlab 和 R 翻译的代码与 @hrbrmstr 的 github hilbert 包进行比较

从 hrbrmstr github hilbert 包中获取 hilbert.cpp 文件

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

将mat_Rcpp 输出与mat_R 和matlab 输出进行比较。它与它们不匹配，因此该包中可能存在错误或提供的matlab代码存在问题。

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 使用 hrbrmstr 的 hilbert 包的 matlab 到 R 翻译代码

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) 为什么要列出退货清单？ (b) 将其与 github.com/hrbrmstr/hilbert 进行比较会很有趣

忘了补充一下，对于这样的操作，这是一个很好的翻译示例。对于希望用其他算法做类似事情的其他人来说，这应该是一个很好的参考。我还想在所有裸数字的末尾添加一个L，以确保它们保持为整数，以实现最快的数学运算。

@Sathish 是的，我当然非常感激。一件事：在您进行编辑后，尤其是完整性检查后，我收到错误消息：all(sapply(as.list(match.call()), is.integer)) 不正确。这是我正在做的事情： ord

@Sathish 我已将 i 和 ord 都设置为整数（使用 as.integer()），并且我在顶部包含了 set.seed(1L)，但仍然出现错误。另一个问题：请注意，矩阵与 NA 接壤。 d2xy() 不应该为 2^ord-by-2^ord 矩阵的所有单元格输出坐标吗？

@hrbrmstr 我认为您的 hilbert 包中可能存在错误。请查看我编辑的答案并确认。谢谢