如何将按位函数从 matlab/C 翻译成 R?特例:希尔伯特曲线算法
How to translate bitwise functions from matlab/C to R? Particular case: Hilbert curve algorithm
我正在尝试将用 matlab 编写的脚本转换为 R。该脚本根据希尔伯特曲线将一维坐标映射到二维坐标。
脚本中有一行我不确定如何翻译成 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 逻辑门在原则上是什么。
如果有人翻译了整个剧本,我不会抱怨。
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
输出。与他们不匹配,所以可能是这个包有bug,或者提供的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 matlab 到 R 使用 hrbrmstr 的 hilbert 包翻译的代码
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
我正在尝试将用 matlab 编写的脚本转换为 R。该脚本根据希尔伯特曲线将一维坐标映射到二维坐标。
脚本中有一行我不确定如何翻译成 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 逻辑门在原则上是什么。
如果有人翻译了整个剧本,我不会抱怨。
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
输出。与他们不匹配,所以可能是这个包有bug,或者提供的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 matlab 到 R 使用 hrbrmstr 的 hilbert 包翻译的代码
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