将C代码转换为R代码 : parsing to change a C function in R (pow(a, b) to a^b)

Question

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

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

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

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

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

按照早期帖子中关于将 Mathematica 代码转换为 R 代码的建议 (Convert Mathematica equations into R code) ，我知道合理的做法是将 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))"))

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

Answer 1

这里有多个问题:

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

但是，根据您的要求(如果我正确理解您想要实现的目标)和教育目的；这是一个示例，其中我们将使用 R 转换您的“伪 C”字符串并创建一个内联的 cfunction()

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

---

假设以下等式:

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

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

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

创建一个“要重新编码的函数”的命名 vector (以及一个没有它们和数字的子集)

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