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指出的那样,这是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处理等.

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