在R中快速限定数据

Question

假设我有一个向量vec，它很长（从1E8个条目开始），并且想将其限制在范围[a,b]内。 我可以编写vec[vec < a] = a和vec[vec > b] = b，但这需要对数据进行两次扫描，并为临时指示器向量（约800MB，两倍）分配大量RAM。 两次扫描会浪费时间，因为如果我们将数据从主存储器复制到本地缓存一次，就可以更好地处理它（调用主存储器或缓存未命中都不好）。还有谁知道使用多个线程可以改善多少，但让我们别太贪心。 :)

是否有某个基于R的良好实现或一些我忽略的包，或者这是Rcpp（或我老朋友 data.table ）的工作？

- Iterator

Answer 1

一个朴素的 C 解决方案是

library(inline)

fun4 <-
    cfunction(c(x="numeric", a="numeric", b="numeric"), body4,
              language="C")
body4 <- "
    R_len_t len = Rf_length(x);
    SEXP result = Rf_allocVector(REALSXP, len);
    const double aa = REAL(a)[0], bb = REAL(b)[0], *xp = REAL(x);
    double *rp = REAL(result);

    for (int i = 0; i < len; ++i)
        if (xp[i] < aa)
            rp[i] = aa;
        else if (xp[i] > bb)
            rp[i] = bb;
        else
            rp[i] = xp[i];

    return result;
"
fun4 <-
    cfunction(c(x="numeric", a="numeric", b="numeric"), body4,
              language="C")

如 Dirk 所指出的，使用 CFLAGS = -fopenmp在 ~/.R/Makevars 中就可以实现简单的并行版本，前提是平台和编译器支持 openmp。

body5 <- "
    R_len_t len = Rf_length(x);
    const double aa = REAL(a)[0], bb = REAL(b)[0], *xp = REAL(x);
    SEXP result = Rf_allocVector(REALSXP, len);
    double *rp = REAL(result);

#pragma omp parallel for
    for (int i = 0; i < len; ++i)
        if (xp[i] < aa)
            rp[i] = aa;
        else if (xp[i] > bb)
            rp[i] = bb;
        else
            rp[i] = xp[i];

    return result;
"
fun5 <-
    cfunction(c(x="numeric", a="numeric", b="numeric"), body5,
              language="C")

以及基准测试

> z <- runif(1e7)
> benchmark(fun1(z,0.25,0.75), fun4(z, .25, .75), fun5(z, .25, .75),
+           replications=10)
                 test replications elapsed  relative user.self sys.self
1 fun1(z, 0.25, 0.75)           10   9.087 14.609325     8.335    0.739
2 fun4(z, 0.25, 0.75)           10   1.505  2.419614     1.305    0.198
3 fun5(z, 0.25, 0.75)           10   0.622  1.000000     2.156    0.320
  user.child sys.child
1          0         0
2          0         0
3          0         0
> identical(res1 <- fun1(z,0.25,0.75), fun4(z,0.25,0.75))
[1] TRUE
> identical(res1, fun5(z, 0.25, 0.75))
[1] TRUE

在我的四核笔记本电脑上运行。假设输入为数字，没有错误检查和NA处理等。

- Martin Morgan

---

+1 我希望在核心R中有这个函数，类似于 clamp(x, low, high)...好吧，总是可以期望的，对吧 ;-) - Tommy

+1 for OpenMP，但我认为你需要修改PKG_CFLAGS等来获取-fopenmp。或者你在其他地方做了这样的修改，例如在~/.R/Makevars中吗？ - Dirk Eddelbuettel

@DirkEddelbuettel R的configure.ac检测到OpenMP；-fopenmp在R_HOME/etc/Makeconf中设置。 - Martin Morgan

不在我的机器上，我从你的示例中得到了“warning: ignoring #pragma omp parallel [-Wunknown-pragmas]”警告。尽管我在/etc/R/Makeconf中有-fopenmp（它是指向R_HOME下面的位置的符号链接）。 - Dirk Eddelbuettel

@DirkEddelbuettel 是的，你说得对，~/.R/Makevars 包含 CFLAGS = -fopenmp。 - Martin Morgan

Answer 2

首先，你的解决方案与pmin/pmax解决方案几乎没有什么区别（我使用n=1e7来尝试而不是n=1e8，因为我不太耐烦）——实际上，pmin/pmax稍微慢一些。

fun1 <- function(x,a,b) {x[x<a] <- a; x[x>b] <- b; x}
fun2 <- function(x,a,b) pmin(pmax(x,a),b)
library(rbenchmark)
z <- runif(1e7)

benchmark(fun1(z,0.25,0.75),fun2(z,0.25,0.75),rep=50)

                 test replications elapsed relative user.self sys.self
1 fun1(z, 0.25, 0.75)           10  21.607  1.00000     6.556   15.001
2 fun2(z, 0.25, 0.75)           10  23.336  1.08002     5.656   17.605

- Ben Bolker

---

有趣。我原以为那会更快，但看来没有这样的运气。 - Iterator

在我的R版本2.15.0 Patched (2012-05-01 r59304)平台上，x86_64-unknown-linux-gnu（64位）使用CFLAGS=-O0编译，fun2比fun1快约20％，而hack .Internal(pmin(FALSE, x, a))等则比fun1快约30％。 - Martin Morgan