将 C 代码转换为 R 代码:解析以更改 R 中的 C 函数(pow(a,b) 到 a^b)

Converting C code into R code: parsing to change a C function in R (pow(a,b) to a^b)

我正在使用 Mathematica 将方程生成为 C 代码(使用 CForm[]),以便将方程导出为字符串并在 R 中使用它。

例如,作为字符串导入到 R 中的 CForm[] 输出如下所示:

"Tau * Power(Omega * (-(R * Gamma) + R),(Tau + R))"

我的问题是如何最好地将上面的 C 代码转换成这样的 R 表达式:

Tau * (Omega * (-(R * Gamma) + R ))^(Tau + R)

根据之前 post 关于将 Mathematica 代码转换为 R 代码 的建议,我知道合理的做法是将 Power() 重新定义为函数,即,:

Power <- function(a,b) {a^b}

但是,通过一系列测试,我发现评估以下形式的表达式:

eval(parse(text="Tau * (Omega * (-(R * Gamma) + R ))^(Tau + R)"))

比将 Power() 定义为函数并评估以下内容的替代方法快得多(在我的 mac 上大约快 4 倍):

eval(parse(text="Tau * Power(Omega * (-(R * Gamma) + R),(Tau + R))"))

这似乎是一个复杂的模式匹配问题,但我找不到任何解决方案。我很感激任何建议。

这里有多个问题:

  1. 您的方程式不是标准 C 代码。来自 Mathematica 的 CForm[] 没有将您的代码转换为正确的 C 语法。也许你可以按照 this answer 并使用 SymbolicC 来解决这部分
  2. 你的问题更多是关于从语言 A 到语言 B 的解析。正如@Olaf 在评论中提到的:你最好使用真正的 C 函数并从 R 调用它或转换手动,具体取决于您执行此操作的频率

但是,根据您的要求(如果我正确理解您想要实现的目标)并且出于教育目的;这是一个示例,其中我们将使用 R 转换您的 "pseudo-C" 字符串并创建一个内联的 cfunction()

注意:这绝不是为了优雅或实用,但总体思路应该希望能帮助您入门


假设以下等式:

v1 <- "4 * Power(Omega * (-(R * Gamma) + R),(Tau + R))"

从原始字符串中提取所有变量和函数

n1 <- stringi::stri_extract_all_words(v1)[[1]]

创建 "functions to recode" 的命名向量(以及没有它们和数字的子集)

newFunc <- c("Power" = "pow")   
n2 <- setdiff(n1, names(newFunc))
n3 <- n2[is.na(as.numeric(n2))]

构建替换列表以供 gsubfn()。为了这个例子,我们用新函数替换旧函数并用 asReal() 包裹变量

toreplace <- setNames(
  as.list(c(newFunc, paste0("asReal(", n3, ")"))), 
  c(names(newFunc), n3)
)

v2 <- gsubfn::gsubfn(paste(names(toreplace), collapse = "|"), toreplace, v1)

然后您可以将这个新字符串传递给 cfunction() 以在 R

中执行
#install.packages("inline")
library(inline)
foo <- cfunction(
  sig  = setNames(rep("integer", length(n3)), n3), 
  body = paste0(
    "SEXP result = PROTECT(allocVector(REALSXP, 1));
     REAL(result)[0] = ", v2, "; 
     UNPROTECT(1);
     return result;"
  )
)

这应该比使用 eval(parse("..."))^ 或定义 Power() 函数

更快
Tau = 21; Omega = 22; R = 42; Gamma = 34
Power <- function(x,y) {x^y}

microbenchmark::microbenchmark(
  C  = foo(Omega, R, Gamma, Tau),
  R1 = eval(parse(text="4 * ((Omega * (-(R * Gamma) + R ))^(Tau + R))")),
  R2 = eval(parse(text="4 * Power(Omega * (-(R * Gamma) + R),(Tau + R))")), 
  times = 10L
)

#Unit: microseconds
# expr     min      lq     mean   median      uq      max neval
#    C   1.233   2.194   5.9555   2.9955   3.302   34.194    10
#   R1 190.012 202.781 230.5187 218.1035 243.891  337.209    10
#   R2 189.162 191.798 374.5778 207.6875 225.078 1868.746    10