将 C 代码转换为 R 代码：解析以更改 R 中的 C 函数（pow(a,b) 到 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 在 cmets 中所提到的：您可能会更好地使用真正的 C 函数并从 R 调用它或手动转换它，具体取决于关于你这样做的频率

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

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

---

假设以下等式：

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

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

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

创建一个包含“要重新编码的函数”的命名向量（以及一个没有它们且没有数字的子集）

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