R中数据的快速边界

Question

假设我有一个向量 vec，它很长（从 1E8 个条目开始）并且想将其绑定到 [a,b] 范围内。我当然可以编写vec[vec < a] = a 和vec[vec > b] = b，但这需要两次传递数据和为临时指标向量分配大量 RAM（~800MB，两次）。两者都通过了刻录时间，因为如果我们只将数据从主存复制到本地缓存一次，我们可以做得更好（对主存的调用很糟糕，缓存未命中也是如此）。谁知道多线程可以改善多少，但我们不要贪心。 :)

在基本 R 或某些我忽略的包中是否有很好的实现，或者这是 Rcpp（或我的老朋友 data.table）的工作？

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 指出的，这是在 ~/.R/Makevars 中使用 CFLAGS = -fopenmp，并且在支持 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 处理等。

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